Louis De Lary De Latour
email: l.delarydelatour@brgm.fr

Annick Loschetter

Olivier Bouc

Jeremy Rohmer

Curtis M Oldenburg

L De Lary

Assessing health impacts of CO2 leakage from a geological storage site into buildings: Role of attenuation in the unsaturated zone and building foundation

Keywords: Geological storage, Leakage, CO 2, Health, Unsaturated zone, Indoor exposure. f) Main

Introduction

Carbon dioxide capture and storage (CCS) is now widely acknowledged as one of the promising mitigation options for reducing atmospheric greenhouse gas emissions. It provides a way to avoid releasing CO 2 in the atmosphere by capturing CO 2 from industrial facilities, transporting it to a storage site and injecting it into a suitable underground geological formation, typically at depth greater than 800 m. The purpose of geological storage is permanent containment of CO 2 (IPPC, 2005). As any industrial activity, geological storage may lead to potential risks on human health and environment. To be viable, a CO 2 storage project must be environmentally safe and socially sustainable, thus, it must be demonstrated that no significant leak will occur [START_REF] West | Issue Profile: Environmental Issues and the Geological Storage of CO 2[END_REF]. The European Directive 2009/31/EC which establishes a legal framework for the deployment of the geological storage of CO 2 , underlines the fact that "a geological formation shall only be selected as a storage site, [… ] if no significant environmental or health risks exist" (European Commission, 2009). To ensure this goal, the Directive requires the consideration of numerous criteria for the characterization and assessment of the storage complex. In particular, a risk assessment that considers the effects of leakage on environment and human health must be carried out. A risk analysis should take into account the consequences (i.e. the impacts) and the probabilities of occurrence (i.e. the likelihood) of unexpected events.

Among the unexpected events, the migration from the reservoir formation to the surface by a preferential permeable pathway (e.g. a faulted zone or an improperly abandoned well) has been determined to be one of the most critical risk scenarios because the flow rate may be high in such situations [START_REF] Lewicki | Natural and industrial analogues for leakage of CO 2 from storage reservoirs: identification of features, events, and processes and lessons learned[END_REF][START_REF] Oldenburg | Certification framework based on effective trapping for geologic carbon sequestration[END_REF]. To date, most research has focused on the assessment of the CO 2 flow rate that could migrate away from the reservoir (e.g. [START_REF] Nordbotten | Model for CO 2 leakage including multiple geological layers and multiple leaky wells[END_REF]Pruess, 2011a) rather than on the impacts of such a flow [START_REF] Price | The consequences of failure should be considered in siting geologic carbon sequestration projects[END_REF]. However, the effects of CO 2 depend much more on the level of exposure (defined as a quantity of substance put in contact with the body barriers) than on the total quantity of CO 2 released [START_REF] Hepple | Carbon Dioxide Capture for Storage in Deep Geologic Formations -Results from the CO 2 Capture Project[END_REF].

The objective of the present paper is twofold: 1) quantify the impacts on humans of CO 2 leakage from a geological storage site into a building in a more realistic way than what has been done so far and, 2) make the results usable within an integrated risk assessment framework. It is worth mentioning that we focus on impacts and not on risks, since we do not quantify the probability of occurrence (which is likely to be very site-specific).

In the first part of this paper, we investigate the potential health effects of CO 2 (section 2). Then, current approaches to characterize hazardous components migration in the unsaturated zone (i.e., the zone between the ground level and the water table) and intrusion into a building are reviewed in different contexts and in the CCS domain (section 3). After identifying limits of current available approaches, we propose a conceptual model to assess CO 2 indoor concentration and associated hazards resulting from a leakage due to a transmissive faulted zone which emerges at some meters depth under a building (section 4). We first investigate the migration of CO 2 in the unsaturated zone and the foundation through numerical simulation with sensitivity analysis (section 5). In section 6, we discuss these results more specifically regarding the relevance of the considered leakage scenario (likelihood and magnitude) and the assumptions of the modeling procedure. On this basis, we finally show how these results could be integrated in a global, efficient and easy-to-use tool to support decision-making of stakeholders and authorities reviewing risk assessment procedures for CO 2 storage projects.

Potential health effects of CO 2

Carbon dioxide is a naturally occurring gas which exists in ambient air in low proportion (about 0.038 % in volume or 380 ppm). It is harmless and even essential to life at low concentration, but it can become hazardous or even lethal at high concentration. Furthermore, carbon dioxide is a color-free gas that cannot be detected by taste or smell below concentration of 10 to 20 % (NIOSH, 1976). The effects of CO 2 on healthy humans are quite well known. Exposure to CO 2 can be hazardous by two mechanisms: (1) because it reduces the oxygen content of ambient air and thus can lead to asphyxia; (2) because of its own direct toxicity. From [START_REF] Hepple | Carbon Dioxide Capture for Storage in Deep Geologic Formations -Results from the CO 2 Capture Project[END_REF], direct toxicity of high levels of CO 2 occurs before the effects of oxygen lack. Carbon dioxide behaves as a threshold substance. This means that the effects of its inhalation are entirely defined by a concentration and duration of exposure. From 1 to 3 % CO 2 leads to a decrease of blood pH and a stimulation of respiration. Above 3 %, the breathing rate increases exponentially, hearing loss or visual disturbance may occur. When CO 2 concentration exceeds 5 %, symptoms are more severe and unconsciousness can occur after one to several minutes at concentrations reaching 10-15 %. At level reaching 30 % of CO 2 death occurs within a few minutes. In France, public authorities have prescribed CO 2 acute toxicity thresholds to be used for industrial hazard studies (Table 1). The "15 minutes Short-Term Exposure Limit Level" of 3 % and the "Immediately Dangerous to Life and Health level" of 4 % are used in the United States [START_REF] Hepple | Carbon Dioxide Capture for Storage in Deep Geologic Formations -Results from the CO 2 Capture Project[END_REF]. Robinson (2010) considers a cautious "action level" (i.e. level above which measures should be taken) of 0.5 % in building. First lethal effects 10%

Significant lethal effects 20%

On natural sites where strong degassing of deep-origin CO 2 occurs, it has been noted that the CO 2 level in soil could be very high (up to 95 % or even more), while in the atmosphere at an elevation of 1 or 2 meters above ground level the CO 2 content is often only twice to thrice the usual concentration [START_REF] Carapezza | Gas hazard assessment in a densely inhabited area of Colli Albani Volcano (Cava dei Selci, Roma)[END_REF][START_REF] Farrar | Magmatic carbon dioxide emissions at Mammoth Mountain, California[END_REF]. Thus, CO 2 is able to dilute rapidly in the atmosphere in open terrain with adequate ventilation. However, because CO 2 is 50 % denser than air (in normal temperature and pressure conditions), it tends to migrate downwards and accumulate in depressions or in confined space potentially creating much higher concentrations [START_REF] Chow | Modeling the effects of topography and wind on atmospheric dispersion of CO 2 surface leakage at geologic carbon sequestration sites[END_REF][START_REF] Benson | Underground geological storage[END_REF][START_REF] Bogen | Source terms for CO 2 risk modeling and GIS/simulation based tools for risk characterization[END_REF]. At natural strong degassing sites, human fatalities due to CO 2 have been recorded under particular situations [START_REF] Carapezza | Gas hazard assessment in a densely inhabited area of Colli Albani Volcano (Cava dei Selci, Roma)[END_REF][START_REF] Holloway | Natural emissions of CO 2 from the geosphere and their bearing on the geological storage of carbon dioxide[END_REF][START_REF] Rogie | Dynamics of carbon dioxide emission at Mammoth Mountain, California[END_REF]. [START_REF] Annunziatellis | Short-and long-term gas hazard: the release of toxic gases in the Alban Hills volcanic area (central Italy)[END_REF] have shown that a danger may exist for exposure in the basement of houses situated near localized sites where strong natural degassing occurs under specific geological conditions. Human fatalities due to CO 2 have also been recorded in other situations such as work in poorly ventilated space [START_REF] Louis | Intoxication par inhalation de dioxyde de carbone[END_REF].

3 Review of current approaches to assess indoor exposure to hazardous compounds.

Analogies with other domains

In order to investigate the issue of CO 2 intrusion into buildings, we took interest in the literature related to other hazardous compounds for which more data are currently available. Some similarities between CO 2 and volatile organic compounds (VOC) have been established, especially regarding the transport in the unsaturated zone (e.g. [START_REF] Zhang | Unsaturated Zone Remediation of Carbon Dioxide Leakage from Geologic Carbon Dioxide Sequestration Sites[END_REF]. The main difference comes from the consequences on environment and human health, the VOCs being potentially harmful at low concentrations. For these compounds, risk assessment methodologies and remediation measures have been available for more than 20 years [START_REF] Epa | Analysis of Selected Enhancements for Soil Vapor Extraction[END_REF][START_REF] Khan | An overview and analysis of site remediation technologies[END_REF]. Two main entry mechanisms from underlying soil to indoor air are identified: diffusion or advection through the structural elements or through openings in the structural elements. Diffusion is the result of a concentration gradient between the underlying soil and the indoor air. Advection is due to a difference of pressure between the soil and the interior of the building. Pressure differentials may be caused by numerous mechanisms such as thermal differences, wind, barometric pressure changes or building ventilation [START_REF] Patterson | Quantification of Vapor Intrusion Pathways into a Slabon-Ground Building under Varying Environmental Conditions[END_REF]. Field study by [START_REF] Patterson | Quantification of Vapor Intrusion Pathways into a Slabon-Ground Building under Varying Environmental Conditions[END_REF] show that the molecular diffusion of VOC vapors through the concrete of a slab-on-grade is the dominant intrusion pathway under normal meteorological conditions. However, when the interior of the building is under artificially reduced pressure, their results suggest that the advective vapors intrusion through gaps or cracks in the slab becomes dominant. The authors underline that these results are highly dependent on site specific conditions.

Numerous analytical or numerical models have been developed to quantitatively assess the intrusion of radon1 or other contaminants into buildings (e.g. [START_REF] Améon | Experimental and theoretical study of radon levels in a house[END_REF][START_REF] Andersen | Numerical modelling of radon-222 entry into houses: an outline of techniques and results[END_REF][START_REF] Johnson | Heuristic model for predicting the intrusion rate of contaminant vapors into buildings[END_REF]. However, these models are not directly applicable in the case of CO 2 leakage because they use as input data a concentration of contaminant in soil; whereas current approaches for CO 2 calculate a mass flow rate through a flow transport model from the reservoir to a shallow aquifer formation or to the surface (e.g. [START_REF] Chang | Simplified Model for CO 2 Leakage and its Attenuation due to Geological Structures[END_REF][START_REF] Pruess | Leakage of CO 2 from geologic storage: Role of secondary accumulation at shallow depth[END_REF].

Analogies may also be found with CO 2 pollution cases: 1) in case of natural degassing of deep origin CO 2 ; 2) in buildings situated near reclaimed or abandoned coal mines [START_REF] Robinson | Occurrence and attempted mitigation of carbon dioxide in a home constructed on reclaimed coal-mine spoil[END_REF]. However, only few investigations have been made about the characteristics and mechanisms of CO 2 intrusion in this context.

Critical review in the CCS domain

Two main approaches have been suggested to assess the exposure to CO 2 :

1. In case of unexpected migration of CO 2 from reservoir into the atmosphere, numerical dispersion calculations have been suggested in order to assess the exposure level of vulnerable entities in open areas. Simulations of plume dispersion have been performed for point sources of CO 2 such as wells (e.g. [START_REF] Aines | Quantifying the potential exposure hazard due to energetic releases of CO 2 from a failed sequestration well[END_REF][START_REF] Bogen | Source terms for CO 2 risk modeling and GIS/simulation based tools for risk characterization[END_REF][START_REF] Chow | Modeling the effects of topography and wind on atmospheric dispersion of CO 2 surface leakage at geologic carbon sequestration sites[END_REF] and for more dispersed source zones [START_REF] Oldenburg | Coupled Vadose Zone and Atmospheric Surface-Layer Transport of Carbon Dioxide from Geologic Carbon Sequestration Sites[END_REF]. Several difficulties have been underlined by these studies, especially about taking into consideration the risks of CO 2 accumulation in low-lying areas, the characteristics of the CO 2 source, and the extent of the reservoir footprint.

2. In order to quantitatively evaluate the impact on human health at the Weyburn site, [START_REF] Stenhouse | Assessing environmental impacts from geological CO 2 storage[END_REF] proposed a conceptual model of exposure to CO 2 into a dwelling resulting from a leakage through an abandoned well. This model converts a flow of CO 2 into an indoor concentration in a very simple way and with few input data (leakage flow rate, volume of dwelling, building ventilation rate) in comparison to the numerical dispersion approach. The calculation of the indoor CO 2 level into the dwelling is based on a mass balance at steady-state conditions. However, this model is overly conservative because it assumes that the CO 2 migrates from the wellhead to the dwelling without any attenuation or resistance during the transport via the upper soil layers or the foundation of the building. It is equivalent to the unlikely situation where a permeable pathway directly emerges inside the building. Furthermore, simulations by [START_REF] Oldenburg | On leakage and seepage from Geologic Carbon Sequestration Sites: Unsaturated zone attenuation[END_REF] showed that the unsaturated zone has a potential, although limited, to attenuate a CO 2 leakage due to different processes such as solubility trapping, gas ponding due to CO 2 density, or dilution in soil gas.

A model considering the transport of CO 2 leaking from a wellhead to a building was suggested by [START_REF] Duguid | Geologic CO 2 sequestration in abandoned oil and gas fields and human health risk assessment[END_REF]. It assumes that the only mechanism responsible for the transport of CO 2 in the unsaturated zone is diffusion. However, simulations carried out by [START_REF] Oldenburg | Modeling Gas Transport in the Shallow Subsurface during the ZERT CO 2 release test[END_REF] have shown that both advection and diffusion control CO 2 transport in the unsaturated zone. Thus the model suggested by [START_REF] Duguid | Geologic CO 2 sequestration in abandoned oil and gas fields and human health risk assessment[END_REF] is likely to underestimate the CO 2 intrusion into buildings.

In the present paper, we chose to focus on the second type of approach to assess the exposure in an enclosed (low ventilated) space, CO 2 being more likely to be harmful in such a situation [START_REF] Oldenburg | Certification framework based on effective trapping for geologic carbon sequestration[END_REF].

Methods

Proposed general conceptual model

The general conceptual model is presented in Figure 1. We suppose that, due to an insufficient site characterization during site selection or an unexpected behavior of the storage complex, the CO 2 is able to migrate upward and to reach the surface by a leakage pathway. We consider that the leakage pathway is a permeable faulted zone situated under an inhabited area. In sedimentary basins, such as the Paris Basin which could be possible candidate for deep saline aquifers storage [START_REF] Grataloup | A site selection methodology for CO2 underground storage in deep saline aquifers: case of the Paris Basin[END_REF], the parent rock is covered by a layer of regolith and soil. Thus, we consider that the leakage pathway does not emerge directly in the house but at some meters below the soil surface. In a zone having an old history of hydrocarbon exploration and exploitation, the leakage pathway may also be an improperly abandoned or an orphan well which has not been mapped during site characterization [START_REF] Oldenburg | Certification framework based on effective trapping for geologic carbon sequestration[END_REF]. The migration from the reservoir towards the surface could be much more complicated than suggested in Figure 1, because of various phenomena such as secondary accumulation or migration in permeable formations situated on the leakage pathway [START_REF] Chang | Simplified Model for CO 2 Leakage and its Attenuation due to Geological Structures[END_REF][START_REF] Nordbotten | Semianalytical solution for CO 2 leakage through an abandoned well[END_REF][START_REF] Pruess | Leakage of CO 2 from geologic storage: Role of secondary accumulation at shallow depth[END_REF]. Furthermore, under the soil layer, the parent rocks are often highly fractured, and thus the migration near the surface is probably more complex than suggested.

Figure 1 -Conceptual model for CO 2 migration from a reservoir and point-source CO 2 leak in the unsaturated zone under a building.

We focus on what happens when the leaking flow reaches the top bound of a transmissive faulted zone situated some meters under a building. To reach the indoor air of the building, the CO 2 has to migrate in the unsaturated zone and through the foundation. During this migration, several attenuation processes may occur. We propose to study these mechanisms and the associated consequences on human health with the conceptual model presented in the circle of Figure 1. The CO 2 source in the unsaturated zone is situated at some meters depth under a building. The building is a slab-on-grade house whose foundation is composed of a layer of gravel and a layer of concrete. This sort of building is more vulnerable to geological gas intrusion than buildings with crawl-space. The interior of the building is modeled as a global volume with a known ventilation rate. This volume can be considered either as a basement or as a ground-level story, depending on the kind of construction. Note that the CO 2 entry in the unsaturated zone is conceptualized as a point source which is located directly beneath the center of the building: we chose this geometrical hypothesis in order to be in a conservative configuration regarding the CO 2 intrusion. To perform this study, two specific models were successively used: 1) a model of migration of CO 2 in the unsaturated zone that assesses the intrusion rate into the building; 2) a model of exposure which converts the intrusion rate into the building into a volume concentration.

Model for CO 2 migration and intrusion into the building

We used a numerical simulation approach to simulate the migration of CO 2 through the unsaturated zone and the building foundation considering a large range of configurations.

Simulations were carried out using the multiphase and multi-component simulator TOUGH2 [START_REF] Pruess | TOUGH2 user's guide[END_REF] combined with the research module EOS7CA [START_REF] Oldenburg | Modeling Gas Transport in the Shallow Subsurface during the ZERT CO 2 release test[END_REF], which is designed for near-surface applications where the pseudo component air is present [START_REF] Oldenburg | On leakage and seepage from Geologic Carbon Sequestration Sites: Unsaturated zone attenuation[END_REF].

We used an axisymmetric model for the unsaturated zone and the building foundation, with the symmetry axis in the vertical direction Z passing through the center of the building and the leakage source (Figure 2). The leakage source is at the bottom of the unsaturated zone. Numerical dispersion due to the grid size was found negligible for the simulated time period. A constant atmospheric pressure of 10 5 Pa is kept at the top of the system. The lateral model boundary is located at 250 m from the leakage source to ensure that the CO 2 plume is not affected by this boundary, which is held at constant gas-static pressure corresponding to the initial soil moisture-air gravity and capillary equilibrium. In a conservative approach, a no flow condition is used at the bottom of the model in order to force the CO 2 plume to migrate upwards. This hypothesis may be representative of reality if the underlying stratum is made of a parent rock or an impervious layer. The initial state of modeling corresponds to a steady state capillarity-gravity equilibrium obtained through a preliminary modeling.

The other hypotheses of the model are the following:

-The unsaturated zone is composed of homogeneous and isotropic material with no fractures or preferential pathways. -Temperature is assumed to be constant throughout the model and equals 15 °C.

-The leakage source is situated at the bottom of the unsaturated zone on a 1 m² surface disk with a view to model a point source. We consider that the CO 2 arrives in the unsaturated zone in gaseous state.

Model for CO 2 exposure in the building

CO 2 entering the building accumulates depending on the ventilation rate of the building. We evaluate the indoor concentration with a simple analytical model based on a mass balance for CO 2 which is extensively used in the area of pollutant intrusion in buildings (e.g. [START_REF] Loureiro | Three-dimensional simulation of radon transport into house basements under constant negative pressure[END_REF]UNSCEAR, 1988). Assuming an instantaneous and homogeneous mixing of CO 2 in the building, the CO 2 indoor air level may be expressed as (Eq.1):

λ ρ ρ × × × - - = × × 2 2 ) ( CO out in in CO V C C Q dt dC V (1)
Where C in is the CO 2 level in indoor air (% vol.), C out is the CO 2 level in outdoor air (% vol.), λ is the ventilation rate (s -1 ), V is the volume of the building (m 3 ), ρ CO2 is the density of CO 2 (g/m 3 ), Q is the intrusion rate of CO 2 into the building (g/s).

Solving this equation at steady state provides the order of magnitude for the indoor CO 2 concentration (Eq. 2):

λ ρ × × + = ∞ = 2 ) ( CO out in V Q C t C (2)
In a real condition, the ventilation rate could be subject to important variations within a few hours (UNSCEAR, 1988) due to human activities (opening or closing windows) or meteorological conditions (wind, temperature change, pressure change). Furthermore, the ventilation rate could be contrasted within buildings and within rooms in the same building. The following values (Table 2) are used in this paper [START_REF] Patterson | Quantification of Vapor Intrusion Pathways into a Slabon-Ground Building under Varying Environmental Conditions[END_REF][START_REF] Stenhouse | Assessing environmental impacts from geological CO 2 storage[END_REF]UNSCEAR, 1988).

Base case scenario

The hypothetical reference building is assumed to be circular with a floor surface area of 100 m² and a volume of 250 m 3 . These correspond to standard values for a building [START_REF] Andersen | Numerical modelling of radon-222 entry into houses: an outline of techniques and results[END_REF]UNSCEAR, 1988). The thickness of the slab is 0.20 m. Under the slab there is a 0.20 m highly permeable layer of gravel. The building and the unsaturated zone properties in the base case are shown in Table 3. The relative permeability and the capillarity functions are set to different values depending on the considered material (soil; gravel; concrete). The parameters of the van Genuchten (1980) model for the concrete slab are from [START_REF] Monlouis-Bonnaire | Prediction of the relative permeability to gas flow of cement-based materials[END_REF]. The parameters for van Genuchten relation for the gravel layer are set to arbitrary values in order to attribute a lower liquid saturation to the gravel. The diffusion coefficient of CO 2 in the gaseous phase is set to 1.5x10 -5 m²s -1 [START_REF] Oldenburg | On leakage and seepage from Geologic Carbon Sequestration Sites: Unsaturated zone attenuation[END_REF]. Intrinsic permeability and porosity for both the slab and the gravel layer are based on [START_REF] Andersen | Numerical modelling of radon-222 entry into houses: an outline of techniques and results[END_REF] and correspond to a standard house. The liquid saturation is set to 0.30 for the unsaturated zone in the base case. For the unsaturated zone, the parameters correspond to fine sand. A literature review shows that depending on the situation and the modality of CO 2 uprising from the reservoir, the flow rate which reaches the surface following a leakage from the reservoir could be subject to substantial variations. Simulations of leakage through a faulted zone situated near a hypothetical reservoir performed by Pruess (2011a and2011b) and accounting for self-enhancing and self-limiting effects showed that an order of magnitude for the CO 2 outflow at ground level could be 1 g/s/m². Pruess (2008) used a 1.6 g/s flow rate to study the role of secondary accumulation during CO 2 migration by a fault. [START_REF] Chang | Simplified Model for CO 2 Leakage and its Attenuation due to Geological Structures[END_REF] modeled the upward CO 2 migration from fault and attenuation by surrounding permeable layers and found that the flux rate reaching the top of the leakage pathway could range from nearly zero to 70 g/s/m² depending on the situation. A well leakage model used by [START_REF] Oldenburg | Certification framework based on effective trapping for geologic carbon sequestration[END_REF] for a hypothetical geologic CO 2 storage project leads to a leakage flux rate from 5.10 -4 to about 5 g/s/m² depending of the permeability of the pathway. However, the maximum leakage flow rate from well with open casing from reservoir to the land surface could be much higher [START_REF] Aines | Quantifying the potential exposure hazard due to energetic releases of CO 2 from a failed sequestration well[END_REF]. On natural sites where deep geological CO 2 is degassing, fluxes measured at ground level can vary from background values (usually less than 1.10 -4 to 3.10 -4 g/s/m²) to more than 1 g/s/m² near the vent core [START_REF] Beaubien | The impact of a naturally-occurring CO 2 gas vent on the shallow ecosystem and soil chemistry of a Mediterranean pasture (Latera, Italy)[END_REF]Carrapezza et al., 2003;[START_REF] Farrar | Magmatic carbon dioxide emissions at Mammoth Mountain, California[END_REF]).

In the base case setting we chose a leakage rate in the unsaturated zone of 1 g/s (which corresponds to a 10 -5 /year fractional leakage rate from a 3 million tons reservoir). However we performed a sensitivity analysis to take into account very contrasted situations.

Results

Base case

The simulations show that the concentration of CO 2 in the considered building reaches, after 3 months, a steady value which depends on the ventilation rate (Figure 3). The intrusion rate of CO 2 into the building is about 0.17 g/s in the base case. Comparison with the total leakage rate in the unsaturated zone (1 g/s) shows that only a small proportion of the leaking CO 2 enters the building due to attenuation mechanisms that are discussed below. Hazardous levels from Table 1 are not reached in the base case whatever the ventilation rate. Not considering those attenuation mechanisms leads to an indoor concentration above the French irreversible effect threshold of 5 % for the lowest ventilation rate.

Sensitivity analysis

The CO 2 footprint of storage reservoirs will generally be very large (order of magnitude in the hundreds of km²), and the properties of the unsaturated zone and of buildings in such a wide area are expected to significantly vary. An unforeseen leakage by a fault or an abandoned well may eventually occur anywhere in the footprint of the storage reservoir and thus may occur in very different situations. To be confident that a wide range of situations (including the most deleterious) has been considered, we performed a sensitivity analysis considering a large range of parameters.

From the base case configuration, we performed a single-factor sensitivity analysis based on the OAT technique, which consists in varying one factor at a time to show the influence of this factor [START_REF] Campolongo | Screening methods, in Sensitivity Analysis[END_REF]. The ventilation rate was set to the lowest value (0.13 h - 1 ) (except for the sensitivity analysis on the ventilation rate). The leakage rate value was kept below 20 g/s in order to avoid excessive pressure build-up in the unsaturated zone. The investigated range of permeability for the unsaturated zone corresponds to soil textures varying from gravel to silt. Very wide ranges of permeability and porosity were chosen to model different slab properties. Results are presented on the "Tornado-type" diagram in Figure 4.

Figure 4 -Results of the sensitivity analysis. The values on the figure are the values of the considered parameters for which the extrema are reached. Note that the scale for the intrusion rate is not applicable for the sensitivity analysis on ventilation rate because in this particular case the intrusion rate is the same whatever the situation. The French threshold for irreversible effect due to acute toxicity (5% CO 2 ) is indicated in the figure.

Two main results can be drawn from this sensitivity analysis:

-For the range of parameters used in this sensitivity analysis, the CO 2 indoor concentration varies from nearly zero to about 6 %, which is above the French value for first irreversible effects of 5 %. The Immediately Dangerous to Life or Health level of 4 % is exceeded only for the situations with a low unsaturated zone thickness or a very high flow rate. For all other situations, acute toxicity thresholds are not reached. Therefore, intrusion of leaking CO 2 in low ventilated parts of buildings (basement) may lead, only under very specific conditions, to hazardous situations. -The most influential parameters on the indoor concentration are ranked based on the width of the variation in the following order: 1) the leakage flow rate; 2) the thickness of the unsaturated zone; 3) the permeability of the unsaturated zone; 4) the building ventilation rate. It is worth noting that we found only a limited number of crucial parameters.

Dominant CO 2 transport mechanisms

In our simulations, the transport of CO 2 in the liquid phase of the unsaturated zone is negligible compared to the transport in the gaseous phase. The leakage provokes a pressure build-up in the unsaturated zone ranging in the base case (Figure 5a) from a hundred of Pa at one meter depth under the building to about 4x10 3 Pa at the leakage source zone. With a low permeability (k=10 -14 m²) unsaturated zone (Figure 5b), the pressure build-up is about 4x10 3 Pa whatever the depth under the building. Such variations of pressure are very significant compared to other variations of pressure (due to barometric pumping, thermal effects etc. -cf. discussion in section 6.1) which are known to have strong effects on the intrusion of geological gases. With a high permeability unsaturated zone (k=10 -9 m²), the pressure buildup reaches only a few Pa close to the leakage source zone (not shown). In this situation, the intrusion in the building is very weak and the resulting indoor concentration is only 700 ppm or 0.07 % (Figure 4) due to a very limited advective migration. In the base case the transport is mainly advective near the source of CO 2 where the pressure build-up is the most significant whereas it is essentially diffusive near the soil surface due to a strong concentration gradient and low pressure gradient (Figure 6). These findings are in agreement with the simulations performed on the CO 2 release test at the Montana site [START_REF] Oldenburg | Modeling Gas Transport in the Shallow Subsurface during the ZERT CO 2 release test[END_REF].

Attenuation processes of the leakage rate

In the base case, the flow intruding into the building is significantly attenuated compared to the flow penetrating the unsaturated zone: only about 17 % of the leaking flow (1 g/s) enters the building. We suggest two main mechanisms to explain this attenuation: 1) the horizontal spreading effect of the CO 2 plume in the unsaturated zone, which enables the CO 2 to bypass the foundation of the building; and 2) the barrier effect of the foundation of the building. The horizontal spreading of a CO 2 plume in the unsaturated zone is a phenomenon that has already been underlined by [START_REF] Oldenburg | On leakage and seepage from Geologic Carbon Sequestration Sites: Unsaturated zone attenuation[END_REF]. The barrier effect of foundation has been largely studied in the field of pollutant intrusion into buildings because it can have a strong influence on the component intrusion rate (e.g. [START_REF] Patterson | Quantification of Vapor Intrusion Pathways into a Slabon-Ground Building under Varying Environmental Conditions[END_REF]. We assessed the barrier effect of the foundation through numerical simulations with hypothetical buildings without foundation: it seems to be limited and leads to a less substantial attenuation than the spreading effect. Thus, properties of the foundation are not among the most influential parameters controlling the intrusion of CO 2 (Figure 4).

The spreading effect explains the results of Figure 7: the higher the total flow rate, the larger the plume, and the lower the percentage of the total flow entering the building due to bypassing. Thus, the indoor exposure increases much slower than the leakage rate (Figure 7).

Figure 7 -Influence of the leakage flow rate in the unsaturated zone on the percentage of CO 2 entering the building and on the indoor CO 2 concentration.

Influence of the unsaturated zone properties on CO 2 migration

The permeability of the unsaturated zone also appears to control the intrusion rate. The shape of the plume with contrasted permeability values for the unsaturated zone (assumed homogeneous and isotropic) is shown in Figures 8a, 8b and 8d. In a highly permeable soil, the vertical transport of the plume is limited and the horizontal transport is dominant. Therefore a significant proportion of the flow rate passes around the building and the resulting indoor concentration is low. Figure 4 shows that the indoor concentration increases from 0.3 % to nearly 5 % when the thickness of the unsaturated zone varies from 20 to 1 meter. Again, this difference is due to the spreading of the plume in the unsaturated zone: it becomes less important as its thickness decreases. (cf. Figure 8c).

6 Discussion and perspectives

Discussion concerning the model

The choice of our model aims at obtaining a compromise between high consequences / very low probability cases (e.g. orphan leaking well directly emerging in the building as considered by [START_REF] Stenhouse | Assessing environmental impacts from geological CO 2 storage[END_REF]) and moderate consequences / higher probability (e.g. leakage at hundreds of meters to kilometers form habitations). Thus, this model is less conservative in terms of risk assessment than the one proposed by [START_REF] Stenhouse | Assessing environmental impacts from geological CO 2 storage[END_REF], however, we consider that this model is still conservative for several reasons.

-The position of the building right above the leakage pathway is the most conservative and does not seem very realistic. Additional simulations should be carried out varying the relative position of the leakage pathway and the building; however the model, assumed axisymmetric so far, should be adapted to that purpose. -The building considered is a slab-on-grade building, which tends to reinforce the entry of geological gases. [START_REF] Groves-Kirkby | Domestic radon remediation of U.K. dwellings by sub-slab depressurisation: evidence for baseline contribution from constructional materials[END_REF] show that other situations such as the ventilated sub-slab configuration significantly reduce the intrusion of geological gases into buildings. Furthermore, in modern slab-on-grade foundation, a polyethylene film prevents the intrusion of gases. -In our model the pressure in the building is held constantly equal to atmospheric pressure. However, indoor pressure is often a few Pa lower [START_REF] Andersen | Numerical modelling of radon-222 entry into houses: an outline of techniques and results[END_REF][START_REF] Johnson | Heuristic model for predicting the intrusion rate of contaminant vapors into buildings[END_REF] because of the wind effect and of the "stack effect" due to temperature differences between indoor and outdoor air [START_REF] Améon | Experimental and theoretical study of radon levels in a house[END_REF]. This difference of pressure may lead to an advective gas flow from the soil to the building. However, from [START_REF] Andersen | Numerical modelling of radon-222 entry into houses: an outline of techniques and results[END_REF] and [START_REF] Loureiro | Three-dimensional simulation of radon transport into house basements under constant negative pressure[END_REF] while the permeability of the soil under the foundation is lower than about 10 -11 m², the main mechanism of intrusion remains diffusion, so neglecting in our model indoor depressurization and cracks in foundation may, only for the highest permeability soils, lead to underestimating the proportion of CO 2 intrusion into the building. -From simulation results [START_REF] Oldenburg | On leakage and seepage from Geologic Carbon Sequestration Sites: Unsaturated zone attenuation[END_REF] and on-field experiments [START_REF] Robinson | Occurrence and attempted mitigation of carbon dioxide in a home constructed on reclaimed coal-mine spoil[END_REF], barometric pumping (due to cyclic natural variations of barometric pressure) may cause local maxima both in the CO 2 flux released from the soil and in the resulting exposure concentration in the atmosphere or in buildings. As barometric pumping is not included in our model, the calculated concentration in the building remains constant after a near steady state is reached, whereas in real situation it should fluctuate with time. Thus, indoor CO 2 exposure could be higher than calculated for short periods of time during rapid drops in barometric pressure. However, [START_REF] Oldenburg | On leakage and seepage from Geologic Carbon Sequestration Sites: Unsaturated zone attenuation[END_REF] found that barometric pumping should have negligible effect on the time-average flow of CO 2 which migrates out of the ground. -We performed some additional simulations to show the influence of permeability anisotropy in the vertical direction to take into account hypothetical soil layering.

When the anisotropy factor between radial and vertical permeabilities varies from 1:1 (base case) to 1000:1 the percentage of CO 2 flow entering the building decreases from 17 % to 4 % because of the spreading of the plume in the unsaturated zone.

-Other attenuation processes such as dissolution in a saturated zone or dissolution into the downward infiltrating meteoric water [START_REF] Oldenburg | On leakage and seepage from Geologic Carbon Sequestration Sites: Unsaturated zone attenuation[END_REF] have not been modeled in our conceptual model. In a real situation, such attenuation processes would decrease even more the percentage of CO 2 entering the building and the resulting indoor concentration. However, if rainfall is intense enough, the compression of soil gases by the wetting front may, by contrast, increase temporarily the CO 2 intrusion and the indoor concentration [START_REF] Guo | Rain-induced subsurface airflow and Lisse effect[END_REF][START_REF] Robinson | Occurrence and attempted mitigation of carbon dioxide in a home constructed on reclaimed coal-mine spoil[END_REF].

Moreover, our model is not meant to predict real consequences of a leakage at a CO 2 storage site: the intent is to use it as a first approach to conservatively evaluate what could be the consequences in the case of a leak. As part of an iterative risk management process, risk assessment is used to design appropriate mitigation and monitoring measures. If a leakage pathway is suspected, additional characterization should be performed to properly understand the potential for a leak; the operations should be designed so as to reduce as far as possible the risk of leakage, including programming corrective measures to remediate any detected leak; and a monitoring plan should be implemented to detect any leak before it leads to adverse consequences. We emphasize that a site would not be selected for CO 2 storage if the risk of leakage were found significant. Therefore the consequences computed in our work should not happen; our model is simply one of the tools to be used in this process of assessing the risk related to CO 2 leakage in order to manage it.

Further work is necessary to improve and complete our approach, especially:

-The probability to encounter the modeled scenario leading to hazardous exposure should be assessed. Anyway, in a carefully selected storage site this probability is expected to be very low, as it implies the conjunction of several unlikely circumstances regarding CO 2 leakage and the presence of an inhabited building: (i) the plume has to reach a preferential pathway (faulted zone or abandoned well), (ii) the CO 2 has to migrate along this pathway from depth to the surface in sufficient quantities, (iii) a building has to be built at the top of this conductive pathway (iiii) somebody has to be exposed under specific conditions. -In the proposed approach, we relied on a One-At-a-Time approach for the sensitivity analysis. The future direction of the present work is to estimate quantitatively the importance of each parameter accounting for the possible interactions and correlations between them with more sophisticated tools for a global sensitivity analysis (e.g. [START_REF] Saltelli | Global Sensitivity Analysis: The Primer[END_REF]. -Our modeling approach would benefit from validation by comparison to fields experiments (e.g. monitoring an instrumented building submitted to a CO 2 artificial flow).

Further developments towards a decision-making tool

The scenario considered above involves different physical phenomena which correspond to different spatial zones:

-Firstly, migration in the saturated zone through a fault or a well (not investigated in this paper); -Secondly, migration in the unsaturated zone and the foundation of the building; -Thirdly, accumulation in the building.

In order to use the results presented in this paper in a risk assessment workflow, there is necessity to integrate them in a global tool that compiles various information and knowledge. The need for different process models for the different zones means that it is not possible to carry out a single global simulation from reservoir to ground surface. Therefore we suggest an adequate modular approach (cf. Figure 9) composed of:

-Module 1: CO 2 leakage from the reservoir and migration in the saturated zone through a fault or a well, -Module 2: CO 2 migration in the unsaturated zone and the foundation of the building, -Module 3: CO 2 accumulation and exposure in the building, -Module 4: Exposure to impurities.

Justification of the segmentation

The segmentation in "saturated zone", "unsaturated zone" and "building" may appear obvious given the fact that specific models exist for these different parts. However the relevance of such segmentation with relevant inputs and outputs had to be validated, especially for module 1 and module 2.

Most models concerning the migration of CO 2 plumes in the saturated zone deliver a leakage rate. However, the few meters of soil and the building foundation may behave as a "stopper" and decrease the total leakage rate. This would make the previous segmentation nonrelevant, since there would be a feedback from module 2 on module 1. We carried out an extensive survey based on an analogy between Darcy's law and Ohm's law (e.g. [START_REF] Endo | Analysis of fluid permeation through a particle-packed layer using an electric resistance network as an analogy[END_REF] in order to compute a "resistance parameter". This analytical approach was compared against numerical results and was found to be correct. Using this approach enabled to compare the order of magnitude for resistance in the saturated zone and in the unsaturated zone. Sensitivity analyses were carried out and gave evidence that the resistance due to the unsaturated zone is less, or even greatly less than the resistance in the saturated zone. In a conservative approach, it is then concluded that module 1 and module 2 can indeed be tackled separately and linked together through the leakage rate (output of module 1 and input of module 2).

Module 1: Migration of CO 2 from a geological reservoir through a fault or a well

The development of this module is not the aim of this paper. In order to elaborate a first complete and functional version of the toolbox, we implemented a very simple analytical model based on Darcy's law. A literature review shows a number of models for the migration of CO 2 in the saturated zone that could be used in order to further develop this first module, in particular: [START_REF] Pruess | Numerical Simulation of CO 2 Leakage from a Geologic Disposal Reservoir, Including Transitions from Super-to Sub-Critical Conditions, and Boiling of Liquid CO[END_REF][START_REF] Pruess | Modeling CO 2 leakage scenarios, including transitions between super and sub-critical conditions and phase change between liquid and gaseous CO 2[END_REF], [START_REF] Pruess | Multiphase flow dynamics during CO 2 injection into saline aquifers[END_REF], [START_REF] Chang | Simplified Model for CO 2 Leakage and its Attenuation due to Geological Structures[END_REF], [START_REF] Silin | A Model of Buoyancy-Driven Two-Phase Countercurrent Fluid Flow[END_REF][START_REF] Hayek | Modeling vertical stratification of CO 2 injected into a deep layered aquifer[END_REF], [START_REF] Nordbotten | Semianalytical solution for CO 2 leakage through an abandoned well[END_REF].

Module 2: Migration of CO 2 in the unsaturated zone and through the foundation

This module aims at obtaining the intrusion rate in the building given a point source leakage rate a few meters deep below the building, as described above. The results of the sensitivity analysis may be used in order to obtain an accurate value for a known situation. We suggest that a 60% attenuation of the CO 2 flow rate entering the unsaturated zone constitutes a conservative value for assessing a large area containing a high variety of settings.

Module 3: Exposure in the building

The goal of this module is to convert the intrusion rate in the building into indoor air concentration. We use the analytical model presented above (cf. 4.3). More complex models developed taking into account the indoor circulation between rooms and floors may alternatively be used to develop this module. To broaden the field of application of the tool, it may also be interesting to develop a surface module corresponding to free atmosphere based on atmospheric surface-layer advection and dispersion (cf. [START_REF] Oldenburg | Coupled Vadose Zone and Atmospheric Surface-Layer Transport of Carbon Dioxide from Geologic Carbon Sequestration Sites[END_REF][START_REF] Chow | Modeling the effects of topography and wind on atmospheric dispersion of CO 2 surface leakage at geologic carbon sequestration sites[END_REF][START_REF] Bogen | Source terms for CO 2 risk modeling and GIS/simulation based tools for risk characterization[END_REF].

Module 4: Integration of a realistic CO 2 stream (CO 2 and gaseous impurities)

Leakage from a geological reservoir is likely to be mostly composed of CO 2 . However, the injected CO 2 may contain other components, which we call impurities. Some substances may also be mobilized during the ascent of CO 2 . As shown in Table 4, exposure limits for some gaseous impurities are far lower than for CO 2 . To take into account these impurities in our tool, we propose an independent module, while making strong hypotheses. We suppose that the leaking flow entering into the building has the same composition than the injected CO 2 . This neglects the possible chemical reactions or mobilizations during CO 2 migration in the reservoir and then through upper layers. We also assume that impurities behave as passive gases: due to their low concentrations, their transport is governed by the CO 2 movement. This approach completes our model and enables to obtain orders of magnitude, even though it should be improved in the future and gain more integration with other modules.

Presentation of the global toolbox

The aggregation of these 4 modules into a global toolbox provides an efficient, quick and easy-to-use tool for decision support, which can be used to assess the impact to human health of a leak from a CO 2 storage site, with respect to threshold concentrations. In the current development of the toolbox, shown in Figure 10 with values corresponding to a hypothetical site, the user interacts through an interface implemented in a simple calculation sheet by filling yellow cells concerning the parameters of the leakage pathway and the reservoir (Module 1), the attenuation in the unsaturated zone (Module 2), the building characteristics (Module 3), the composition of the injected gas (Module 4). Then, the conclusions of the impact analysis are immediately generated. In the example Figure 10, effects on heath are possible for H 2 S, because the injected gas contains significant proportion of this highly toxic gas.

Conclusion

We set up a conceptual model to assess indoor exposure to CO 2 resulting from a leakage from a faulted zone or an abandoned well, which is located beneath a building at some meters depth. The aim of this approach is to assess the health consequences on humans in case of leakage. The probability to be exposed to hazardous concentrations is not assessed, but it is expected to be very low in a carefully characterized, assessed and monitored site.

We performed a sensitivity analysis regarding the migration of CO 2 in the unsaturated zone in order to understand the possible attenuating roles of the unsaturated zone and the foundation. For the ranges of parameters investigated, the indoor CO 2 concentration varies from no effect level to hazardous concentrations. A main finding of this study is that hazardous concentrations are reached in a building only under specific conditions: very low indoor ventilation rate, high flow rate and/or building situated very close to a leaking fault or well. The most influential parameters controlling the indoor CO 2 concentration are the leakage flow rate, the thickness of the unsaturated zone, its permeability, and the building ventilation rate.

The CO 2 leakage flow is attenuated mainly due to the horizontal spreading of the plume in the unsaturated zone with a complementary but weaker attenuation due to the barrier effect of the building foundation. The transport of CO 2 is mainly advective close to the leakage source and diffusive close to the surface.

The main benefit of this model is to assess in a more realistic way than previously suggested approaches the exposure in a building following a CO 2 leakage. A site-specific approximation of the indoor exposure could be easily deduced from the most influential parameters. In the general case, a 60% attenuation of the CO 2 flow rate entering the unsaturated zone is considered a conservative value. Nevertheless, further improvements are necessary especially regarding the control of building-related processes on CO 2 intrusion.

We have integrated these results in a global toolbox, so as to be able to assess the consequences on human of a leakage from the reservoir to the surface. The problem has been segmented in three modules corresponding to the different spatial zones: the saturated zone, the unsaturated zone, the building. Their association delivers the CO 2 content in the indoor air. A fourth module roughly estimates the concentrations of impurities, based on a dilution model. The global toolbox provides an efficient, quick and easy-to-use tool for decision support, which can be used to assess the impact to human health of a leak from a CO 2 storage site.

It is worth mentioning that this toolbox deals with impacts and not with risks: if it were to be employed to assess risk, the impacts should be multiplied by the probability of the event.

Benefit could be taken from this computing time efficient tool in a risk assessment process were numerous simulations would be necessary to take into account ranges of uncertainty. Moreover, such a global tool would be useful for decision-makers in a preliminary risk analysis and would emphasize the most critical points that would require more detailed studies. It could for instance be valuable to help authorities and stakeholders auditing risk assessment studies [START_REF] Bouc | Determining safety criteria for CO 2 geological storage[END_REF][START_REF] Oldenburg | Certification framework based on effective trapping for geologic carbon sequestration[END_REF].

Figure 2 -

 2 Figure 2 -Mesh used for the base-case simulation. Here with a 10 m thick unsaturated zone.

Figure 3 -

 3 Figure 3 -Concentration in the building as a function of time in the base case for different ventilation rates.

Figure 5 -

 5 Figure 5 -Pressure increase in the model when a nearly steady state is reached (a) in the base case, (b) with a 10 -14 m² (silt) and 10 m thick unsaturated zone (c) with a 10 -12 m² and 5 m thick unsaturated zone; dP is the pressure variation in Pa.

Figure 6 -

 6 Figure 6 -Advection vs diffusion of the CO 2 as a function of depth at the center of the model in the base case.

Figure 8 -

 8 Figure 8 -Plume of CO 2 in the model when a nearly steady state is reached (a) in the base case, (b) with a 10 -14 m² (silt) and 10 m thick unsaturated zone, (c) with a 10 -12 m² and 5 m thick unsaturated zone, (d) for a 10 -9 m² permeability and 10 m thick unsaturated zone. Xg(CO 2 ) is the mass fraction of CO 2 in the gas phase.

Figure 9 -

 9 Figure 9 -Proposed segmentation in modules

Figure 10 -

 10 Figure 10 -Overview of the global toolbox. Input values must be entered in yellow cells. Figures in rounded frames indicate intermediary results (output of the module below, input of the module above).

Table 1 -Thresholds for CO 2 acute toxicity in France used for hazard studies (Ministry of Ecology, 2007).

 1

	Thresholds	Corresponding CO 2 level (%
		vol.)
	Irreversible effects	5%

Table 2 -Range of building ventilation rates used in our model.

 2

	Ventilation rate (h -1 )
	Basement	0.13
	Building (Low)	0.5
	Building (Average)	1
	Building (High)	2

Table 3 -Properties of the unsaturated zone and the building foundation in the base case setting.

 3

	Material	Thickness	Permeability	Porosity	Capillary pressure,
		(m)	(m²)	(%)	relative permeability
	Concrete slab	0.2	10 -15	20 %	Van Genuchten; Slr=0.25; Sgr=0.01; α=10 -7 ; m=0.5
	Gravel layer	0.2	5x10 -9	40 %	Van Genuchten; Slr=0.1; Sgr=0.01; α=10 -3 ; m=0.4
	Unsaturated zone	10	10 -12	30 %	Van Genuchten; Slr=0.15;Sgr=0.01; α=10 -4 ; m=0.2

Table 4 : Thresholds considered as references for irreversible effects for different gases based on a 30 minutes exposure in France (Baulig et al., 2004; Tissot and Pichard, 2000, 2004a, 2004b)

 4 Consequently, regarding health effects, CO 2 may not be the most critical component. Table5illustrates that a CO 2 flux from precombustion capture (upper value from IPPC, 2005) may contain 165 times less H 2 S than CO 2 ; since the effect threshold is 500 times greater for CO 2 than H 2 S, in a dilution calculation the most penalizing substance would be H 2 S.

	Component CO 2	H 2 S	NO 2	NO	SO 2
	Thresholds 5,00 %	0,010 %	0,005 %	0,010 %	0,010 %

Table 5 : Comparison of ratios concerning maximum exposure concentration and possible composition of injected gas for CO 2 and H 2 S

 5

	CO 2	H 2 S

  Based on these hypotheses we calculate a dilution rate for CO 2 :This dilution rate is then used to obtain the concentration of each gaseous impurity:

	Concentrat	ion	in	the	building	=	Rate injected the Dilution in ion Concentrat	flow
									Ratio 500
									Ratio 165
	Dilution	Rate	=	ion of ion Concentrat Concentrat	CO in 2 CO of 2	building flow injected the the in

Radon is a naturally occurring gas which is present everywhere at the Earth crust. Contrary to CO

radon may lead to adverse impact at very low concentration. Very extensive literature is available about the entry mechanisms of this gas into buildings.


Rémi Blandino
email: remi.blandino@institutoptique.fr

Marco G Genoni
email: m.genoni@imperial.ac.uk

Jean Etesse

Marco Barbieri

Matteo G A Paris

Philippe Grangier

Rosa Tualle-Brouri

Homodyne estimation of Gaussian quantum discord

We address the experimental estimation of Gaussian quantum discord for two-mode squeezed thermal state, and demonstrate a measurement scheme based on a pair of homodyne detectors assisted by Bayesian analysis which provides nearly optimal estimation for small value of discord. Besides, though homodyne detection is not optimal for Gaussian discord, the noise ratio to the ultimate quantum limit, as dictacted by the quantum Cramer-Rao bound, is limited to about 10 dB.

Quantum correlations are central resources for quantum technology. These tight connections empower the advantages shown by the exploitation of quantum coding in applications to cryptography, computation and sensing. While at first entanglement was recognized to be the most peculiar form of quantum correlations, novel concepts have been introduced to capture either more specific aspects, such as quantum steering [1,2], or, to the other end of the spectrum, more general occurrences. Quantum discord represents the most succesful attempt to such a broad picture within the current picture [3,4]. It has recently attracted considerable attention, due to its possible, yet controversial, usefulness as a resource in mixed-state quantum computing [5,6], where entanglement is shown to be exponentially small [7,8]. Furthermore, it has been demonstrated that discord does play a role in the activation of multipartite entanglement [9], entanglement generation by measurement [10], state merging [11], and for complete positivity of evolutions [12,13].

Experimental observation of quantum discord has been undertaken either by direct inspection of the density matrix [14][15][16][17], or by using a witness [18] with no concern about the optimality of the scheme. A key problem is then to find strategies for the quantification of these resources with precision, and to understand its fundamental limit introducing proper Cramér-Rao bounds (CRB) [19][START_REF] Braunstein | [END_REF][21]. This has already been done for the entanglement [22] and optimal estimators have been experimentally proved to attain the quantum limit for different families of qubit states [23]. For the perspective of quantum metrology, this is highly nontrivial, since there exists no observable directly related to quantum discord. A proper estimator is then needed, which might depend on several characteristic parameters of the quantum state. In such a multiparameter problem, finding an optimised detection scheme might be hard, and could demand complex experimental apparata or heavy post-processing of the data.

In this Letter we demonstrate homodyne estimation of Gaussian quantum discord in continuous variable systems [24,25], and compare the achieved level of precision with the classical CRB for homodyne detection, and with the quantum CRB, which sets the ultimate precision allowed by quantum mechanics. Our findings also show how a suitable Bayesian data processing may be employed to improve precision, especially in the estimation of small values of discord.

Our investigation is concerned with an important class of Gaussian states, i.e. the two-mode squeezed thermal states (STS) naturally produced by a non-collinear optical parametric amplifier (OPA). If we introduce the two-mode squeezing operator S 2 (s) = exp s(a † 0 a † 1a 0 a 1 ) , and the thermal state ν

(N )= 1 N +1 n ( N N +1
) n |n n|, we can write the STS as

̺(N s , N t ) = S 2 (s)ν(N t ) ⊗ ν(N t )S 2 (s) † , (1)
and thus can be fully described by the two parameters N s = sinh 2 s and N t , representing, respectively, the effective amount of squeezing photons and thermal photons. In fact, spurious effects such as unwanted amplification, result in a loss of purity of the squeezed state by thermalisation, but do not affect the Gaussian character of the emission, so the form of the density matrix (1) provides a fully general description of the output of a realistic OPA [26].

The study of quantum discord for continuous variable systems has been addressed so far uniquely to Gaussian states and measurement; there exist a general closed formula for discord given an arbitrary covariance matrix [25], but in the following we will use a simplified expression holding for a particular class to which STS belong [24]. We can estimate the discord from N s and N t as we varied the pump power of our OPA [27]. For each power setting, these two parameters are extracted by the outcome of two homodyne detectors, one on each mode, which measure pairs of quadratures {X 0 , X 1 } and {P 0 , P 1 } (Fig. 1). From these, we can evaluate the four linear combinations

Q (1/2) = X 0 ± X 1 √ 2 Q (3/4) = P 0 ± P 1 √ 2 (2)
where Q (1) and Q (4) are squeezed quadratures, while Q (2) and Q (3) are anti-squeezed; in particular M q measurement outcomes are recorded for each one of the four quadratures. The corresponding variances, σ 2 (Q sq ) and σ 2 (Q asq ), that can be obtained from the experimental data, can be rewritten as function of N s and N t as follows

σ 2 (Q sq/asq ) = (1 + 2N s ∓ 2 N s (1 + N s ))(1 + 2N t ) (3)
The expressions obtained can be then inverted to obtain the experimental estimate N inv s and N inv t , along with the relative uncertainties σ 2 (N inv s ) and σ 2 (N inv t ). These values can be used in the expression for discord [27] to calculate its value D inv , and the uncertainty σ 2 (D inv ). The uncertainties on these quantities are then obtained by a Monte Carlo procedure [27]. One can use the same data and refine the estimation by using a Bayesian analysis. As described above, each data sample corresponds to M q = 2 • 10 4 measurement of each of the four quadratures. The total sample, thus correspond to

M T = 4M q homodyne outcomes X = {q (1) 1 , .., q (1) Mq , q (2) 1 , .., q (2)
Mq , q

(3) 1 , .., q (3) Mq , q (4) 1 , .., q (4)
Mq }.

The overall sample probability can be evaluated as

p(X |N s , N t ) = 4 k=1 Mq j=1 p k (q (k) j |N s , N t ) (4)
where the probability of obtaining the outcome q

(k) j by mea- suring the quadrature Q (k) is a Gaussian distribution p k (q (k) j |N s , N t ) = 1 2πσ 2 k exp - (q (k) j ) 2 2σ 2 k . (5)
For squeezed quadratures (k = {1, 4}) we substitute

σ 2 k = σ 2 (Q sq ), while for anti-squeezed quadrature (k = {2, 3}) σ 2 k = σ 2 (Q asq )
. By means of the Bayes theorem, we obtain the a-posteriori probability

p(N s , N t |X ) = 1 N p(X |N s , N t )p 0 (N s )p 0 (N t ), (6)
N = dN s dN t p(X |N s , N t )p 0 (N s )p 0 (N t ). (7)
where the p 0 (N s ) and p 0 (N t ) are the so-called a-priori probability distributions for the two parameters. In our procedure, we use the results of the inversion estimation to construct these a-priori distributions. That is, we consider p 0 (N s ) and p 0 (N t ) as Gaussian functions with respectively, mean values equal to N inv s and N inv t , and variances equal to σ 2 (N inv s ) and σ 2 (N inv t ). Then, we can use the a-posteriori probability distribution evaluated as in Eq. ( 6) to obtain an estimate of the two parameters and of their variances. In formula, (for j = s, t)

N bay j = dN s dN t N j p(N s , N t |X ) (8) σ 2 (N bay j ) = dN s dN t (N j -N bay j ) 2 p(N s , N t |X ) . (9)
By using the formula of the discord for two-mode STS in [27] and by propagating the errors, we then obtain an estimate D bay for the discord, alongs with its variance σ 2 (D bay ).

The value of discord depends on both the squeezing and thermal photons. Consequently, its estimation is inherently a multi-parameter problem, and we have to identify the relevant physical parameters to evaluate the correct CRB. In the multiparametre scenario, the quantum Fisher information (QFI) associated to a vector of parameters λ={λ i } 0≤i≤n is in the form of a matrix H. This sets a lower bound on the covariance σ 2 ij = λ i λ jλ i λ j after M repetitions on the experiment:

σ 2 ij ≥ 1 M H -1 ij (10)
In the specific case of our experiment, we can bound the uncertainty on the discord D of the states we prepare as:

σ 2 (D) ≥ 1 M (H -1
) DD . While our measurement strategy has the advantage of being simple, it is not expected to saturate the quantum CRB; then we also need to compare it to the classical CRB associated to our specific measurement, which is analogously described by a classical Fisher information (FI) matrix F .

In the evaluation of the correct bound, we need a suitable parametrisation of the state, so that in the expression (10) one parameter only actually varies, while the others are kept fixed: this can not be the case for the number of thermal and squeeing photons, as both of them change with the pump power. Therefore, we need to reshape the QFI matrices for different couples of parameters, so to consider those which are more directly connected to the experimental conditions. We start by considering the first couple λ 1 = {N s , N t }; by using the formulas described in [27], we obtain

H (1) = diag (1 + 2N t ) 2 N s (1 + N s )(1 + 2N t + 2N 2 t ) , 1 N t (1 + N t . (11)
As explained above, thermal photons appear because of imperfections in the operation of the OPA and because of loss. When the squeezing is not too low, we can reparametrise our state by taking in consideration the effective squeezing strength r, and a parasite amplification with strength γr [28,29]. The overall homodyne detection can be separately calibrated, obtaining η = 0.62. Thus we can rewrite the matrix (11) in terms of the two unknown physical parameters λ 2 = {r, γ} via the expression

H (2) = B 12 H (1) B T 12
, where B 12 is the transfer matrix for this change of variables [27]. Next, since the physical parameter that changes during our experiment, resulting in the variation of the amount of discord, is the squeezing parameter r (while γ and η can be considered to remain constant), we perform the last change of variable, by considering λ 3 = {D, γ}. Again the QFI matrix can be obtained as

H (3) = B 23 H (2) B T
23 , and the bound on the variance for the quantum discord can be easily evaluated as described in Eq. (10).

We also want to derive the classical CRB for quantum discord, that we obtain if we consider as measurement homodyne detection of squeezed and anti-squeezed quadratures of a twomode squeezed thermal state. Let us start by considering the Fisher information matrix we obtain if we want to estimate the two parameters λ 1 = {N s , N t } by means of homodyne detection on a certain quadrature Q φ . Since the state is a Gaussian state, the conditional probability distribution of measuring a value x, is a Gaussian function, with zero mean, and variance σ 2 (Q φ ). By using the formulas in [27] and evaluaiting some Gaussian integrals, one easily obtains the following formula for the Fisher matrix elements

F µν = 1 2σ 2 (Q φ ) ∂σ 2 (Q φ ) ∂λ µ ∂σ 2 (Q φ ) ∂λ ν (12)
where λ µ = {N s , N t }. If one considers to measure the squeezed or the anti-squeezed quadratures one obtains the following FI matrices:

F sq/asq =   1 2Ns+2N 2 s ∓ 1 √ Ns(1+Ns)(1+2Nt) ∓ 1 √ Ns(1+Ns)(1+2Nt) 2 (1+2Nt) 2  
If we perform a fixed number of measurements, where half of them are done on the squeezed quadratures, and the remaining ones on the anti-squeezed quadratures, the overall FI matrix which will give the CRB for the two parameters λ 1 = {N s , N t } is obtained as

F (1) = F sq + F asq 2 = 1 2Ns+2N 2 s 0 0 2 (1+2Nt) 2 . (13)
To obtain the CRB for the quantum discord, we can proceed as we showed for the quantum CRB, simply replacing the QFI matrices, with the FI ones. The values of the discord obtained using our Bayesian estimation are shown in Fig. 2: the points indicate the experimental data, while the solid line describes the model (1), where the homodyne efficiency η and the relative parasite gain γ are kept to a constant value. Our model is in satisfactory agreement with the data, so we can be confident of that the CRB calculated after the matrix (11) reliably describes the ultimate limit for precision. In the formulae of the Bayes rule ( 6), we need to multiply several probabilities (5), which rapidly give a number hardly manageable by reasonable computing power: this sets a limit to the number of quadrature values one can effectively use in about 800 points.

In order to use larger samples, we have divided our data in N b = 10 2 blocks of 200 points for each of the quadratures (2), calculated the Bayesian estimation of the discord for each block, then considered the average weighted on the associated uncertainties. We notice that the a priori probabilities (5) are calculated from the whole set of data containing M T values: as they intervene in the evaluation for each block, the overall number of resources to be considered is

M = N b • M T .
The comparison between our experimental uncertainties and both the Cramér-Rao limit for our detection (12) and the quantum Cramér-Rao limit is shown in Fig. 3, where we report the quantity K M = M σ 2 (D)/(F -1 ) DD (or the analogue quantity involving the QFI) expressed in dB. K M is the variance of the discord estimator from homodyne data multiplied by the number of resources and divided by the relevant elements of the (quantum) inverse Fisher matrix. For K M equal to unity we have optimal estimation. Solid points refer to Bayesian estimation while empty ones correspond to estimation by inversion. We notice that for low values of discord, the Bayesian technique provides a nearly optimal estimator for the chosen measurement strategy, whereas estimation by inversion is noisier. We also notice that the point corresponding to the lowest value of the discord is slightly below the quantum CRB: this confirms that for low values of the squeezing of the pump, the model we use is not as accurate as in other regimes. For increasing values, the observed variances depart from the optimum by less than an order of magnitude: as Bayesian estimation rapidly converges to optimal, we can attribute this trend to actual variations of the value of the discord in the experiment, becoming more important than statistical fluctuations when the discord increases.

The measurement we have adopted has the considerable advantage of being the simplest experimental option; however, simplicity always comes at a price, and we do not expect it to deliver the best estimator for discord as established by the quantum CRB. In the limit of low discord, we measure a ratio of about 10 dB , which tells us that the price we have to pay is quite reasonable. The departure from the quantum CRB then slightly increases with discord.

In conclusion, we have presented the experimental estimation of Gaussian quantum discord for two-mode squeezed state. Our scheme is based on homodyne detection assisted by Bayesian analysis. Our results are in good agreement with the theorerical model, and this allows us to perform a reliable precision analysis. We found that homodyne estimation shows about 10 dB of added noise compared to the ultimate bound imposed by the quantum Fisher information, with Bayesian analysis that slighlty improves performances for small values of discord. We have also compared our results with the CRB for homodyne detection and found that the estimation is nearly optimal for small values of discord. The usefulness of quantum discord as a resource for quantum technology is a heavily debated topic, and a definitive answer may only come from experiments involving carefully prepared quantum states. Our results contribute to the precise characterization of Gaussian discord and illustrate how a suitable data processing may decrease the uncertainty when optimal detection schemes are not available.

Physical model and evaluation of quantum discord

Here we give explicit expressions of the formulas used in the main text. A two-mode squeezed thermal state (STS) is fully characterized by the two parameters N s = sinh 2 s and N t , representing, respectively, the effective amount of squeezing photons and thermal photons. In our experimental model, these quantities can be obtained as a function of the physical parameters {r, γ, η}, that is

N s = 1 2   -1 + A(r, γ, η) η 2 cosh 4 r cosh 2 (2rγ) + B(r, η) 2 + 2η cosh 2 r(-2η cosh 4 (rγ) sinh 2 r + cosh(2rγ)B(r, η))   N t = 1 2 -1 + (A(r, γ, η) -η cosh 2 (rγ) sinh 2r)(A(r, γ, η) + η cosh 2 (rγ) sinh 2r) (23)
where

A(r, γ, η) = 1 -η + η cosh 2 r cosh 2rγ + η sinh 2 r (24) B(r, η) = 1 -η + η sinh 2 r. (25)
Notice that by varying the pump power, we change the parameter r only, while the noise parameters γ and η stay constant to a very good level of approximation (it should depend on the mode matching only). As a result, both the effective squeezing and thermal photons N s and N t change accordingly.

The covariance matrix of a two-mode STS can be written as

Σ sts = a½ 2 c σ z c σ z a½ 2 (26)
where

a = (1 + 2N t )(1 + 2N s ) (27)
c = 2(1 + 2N t ) N s (N s + 1), (28)
and ½ 2 and σ z are respectively the 2 × 2 identity matrix and the Pauli matrix for the z direction. Following [24], the quantum discord can be thus evaluated, obtaining

expressed as a function of the effective parameters N s and N t . By using Eqs. (23), one can easily obtain the discord as a function of the physical parameters r, γ and η.

Monte Carlo evaluation of uncertainties.

In our experiment, M q measurement outcomes are recorded for each one of the four quadratures. The quadratures Q (1) and Q (4) form a set of 2M q squeezed quadratures measurements, as well as Q (1) and Q (4) form a set of 2M q antisqueezed quadratures measurements. From those two sets of experimental data, we compute the variances σ 2 (Q sq ) and σ 2 (Q asq ). Assuming that those estimated variances follow a Gaussian distribution, the variance of their estimation is given by Var σ 2 (Q sq/asq ) =2σ 4 (Q sq/asq )/(2M q ).

The expressions for sigma 2 (Q sq ) and σ 2 (Q asq ) are then inverted to obtain N s and N t as function of σ 2 (Q sq ) and σ 2 (Q asq ). The mean value and the associated uncertainties are determined by a Monte Carlo simulation of 10 6 experiments. For each experiment, the values σ2 (Q sq ) and σ2 (Q asq ) of the squeezed and anti-squeezed variances are randomly chosen from two Gaussian distributions respectively of mean values σ 2 (Q sq ) and σ 2 (Q asq ), and variances Var σ 2 (Q sq ) and Var σ 2 (Q asq ) .

FIG. 1 .

 1 FIG. 1. Conceptual layout of the discord estimation. Our noncollinear OPA is based on a nonlinear KNbO3 cristal, pumped by a frequency-doubled Ti:Sapph laser (repetition rate 800kHz, wavelength λp=425nm, pulse duration 120fs). This produces a STS with discord D depending on the pump power, i.e. on their average photon number. the two mode are measured by two homodyne detectors DH0 and DH1. The relative phase between the local oscillators is locked by mapping them to different polarisations on the same spatial mode. In this way, we can record blocks of 20000 values of pairs of quadratures for {X0, X1} and the same for {P0, P1}.

FIG. 3 .

 3 FIG.3. The noise ratio KM as a function of discord. Circles and triangles correspond respectively to the quantum and the classical CRB. Solid points refer to Bayesian estimation, while empty ones correspond to estimation by inversion. The uncertainties for the inversion method are estimated using a MonteCarlo procedure on 10 6 points. Notice that the number of resources for Bayesian estimation is M = N b • MT , while for the inversion method M = MT .

  Experimental values of Gaussian quantum discord from homodyne data and Bayesian estimation. The points correspond to the estimated experimental values, while the solid line is the theoretical prediction for η=0.62 and γ=0.73 (the value of γ has been extracted from a best-fit of the points). Uncertainties are within the point size.

	D
	0.10
	0.08
	0.06
	0.04
	0.02
						r
	0.05	0.10	0.15	0.20	0.25	0.30
	FIG. 2.

  D(N s , N t ) = 2N t log(N t ) -2(N t + 1) log(N t + 1) -(N s + N t + 2N s N t ) log(N s + N t + 2N s N t )+ -N t (N t + 1) 1 + N s + N t + 2N s N t log N t (N t + 1) 1 + N s + N t + 2N s N t + (1 + N s + N t + 2N s N t ) log(1 + N s + N t + 2N s N t ) + N s + 2N s N t + (1 + N t ) 2 1 + N s + N t + 2N s N t log N s + 2N s N t + (1 + N t ) 2 1 + N s + N t + 2N s N t

We thank A. Datta, G. Adesso, S. Olivares and M. Paternostro for discussion and comments. We acknowledge support from the EU project COMPAS and the ERANET project HIPERCOM. MGG acknowledges a fellowship support from UK EPSRC (grant EP/I026436/1), MB is supported by the Marie Curie contract PIEF-GA-2009-236345-PROMETEO.

SUPPLEMENTAL MATERIAL

Multiparametric quantum estimation

Here we present the case where the estimation of more than one parameter has to be performed. We define a family of quantum states ̺ λ which depends on a set of N parameters λ = {λ µ }, µ = 1, . . . , N . In this case the geometry of the estimation problem is contained in the QFI matrix, whose elements are defined as

and where we have introduce the Symmetric Logarithmic Derivatives (SLD) L µ corresponding to the parameter λ µ , as the selfadjoint operator that satisfies the equation

In terms of eigenvalues and eigenvectors of ̺ λ , by denoting with ∂ µ the partial derivative respect to λ µ , we have

The QFI matrix here defined provides a lower bound (the quantum Cramér-Rao bound) on the covariance matrix

In the multiparametric case this bound is not in general achievable, on the other hand, the diagonal elements of the inverse Fisher matrix provide achievable bounds for the variances of single parameter estimators, at fixed value of the others

Let us now suppose that we are interested in the estimation of different set of parameters λ = { λν = λν (λ)} which are functions of the previous ones. We then need to reparametrize the family of quantum states in terms of λ. Since ∂ν = µ B µν ∂ µ with B µν = ∂λ µ /∂ λν we have that

and the new QFI matrix simply reads

We consider here the case where we perform a specific indirect measurement in order to infer the values of the parameters λ, given some measurement outcomes X = {x 1 , x 2 , . . . }. The whole measurement process can be described by the conditional probability p(x|λ) of obtaining the value x from the measurement when the parameters have the values λ. Given this object, we can define the Fisher information (FI) matrix whose elements are obtained as

This matrix defines a bound on the covariance matrix γ for the specific measurement we performed. In particular we are interested in the bound for the variance of a single-parameter, at fixed values of the others, which reads

and which in turn is always lower bounded by quantum CRB given in Eq. (18).

Notice that if we have to reparametrize our family of states in terms of different parameters λ, we can use the same formulas shown above for the QFI matrix, obtaining the new FI matrix as F = BF B.

The values of Ñs and Ñt are then computed using those random values. The experimental estimate N inv s and N inv t are finally obtained by taking the mean of the 10 6 values of Ñs and Ñt , whereas their uncertainties σ 2 (N inv s ) and σ 2 (N inv t ) are obtained by computing the variance of the 10 6 values of Ñs and Ñt .

These values can be used in the expression for discord [27] to calculate its value D inv , and the uncertainty σ 2 (D inv ) by using a similar MonteCarlo method.