Jeremy Rohmer
email: j.rohmer@brgm.fr

Daniel Raucoules

On the applicability of Persistent Scatterers interferometry (PSI) analysis for long term CO 2 storage monitoring

Keywords: CO 2 geological storage, Monitoring protocol, InSAR remote sensing, Persistent Scatterers, geomechanical model, Corner Reflectors

Persistent scatterers synthetic aperture radar interferometry (PSI) analysis is viewed as a promising monitoring technique for multiyear detection of ground surface displacements during long term CO 2 injection into deep aquifers (of the order of 25 to 50 years). This technique was successfully utilized at the CO 2 storage site of In-Salah (Algeria).

Nevertheless, this success may not be generalized to future onshore industrial-scale CO 2 storage sites, which are planned in far more complex contexts mostly presenting natural terrains, either agricultural or vegetated areas, i.e. characterized by poor PS density, which might locally drop to zero, hence hampering the deployment of the PSI technique. The objective of the present paper is to discuss the constraints on the PSI application regarding the specificities (magnitude, spatial or temporal scales and evolution) of the deformation phenomenon associated with the CO 2 injection into storage aquifers. This discussion is supported by the estimates of the expected spatio-temporal evolution of surface vertical displacements caused by CO 2 injection in deep aquifers by combining the approximate solutions for pressure build-up and for the associated vertical displacements in an elastic half-1

Introduction

CO 2 capture and geological storage (denoted CCS) is seen as a promising technology in the portfolio of measures required to mitigate the effects of anthropogenic greenhouse gas emissions (IPCC, 2005). Among the available geological targets, deep aquifers are recognized to offer very large potential storage capacity with a broad distribution throughout the world in all sedimentary basins (e.g., IPCC, 2005). Therefore, the present paper will primarily focus on this option. The International Energy Agency (IEA) recently evaluated the contribution of CCS to emissions reductions by 2050 to one fifth in order to achieve the most cost-effective stabilisation of greenhouse gas concentrations (IEA, 2009). On this basis, the IEA concludes that the implementation of the technology should reach 100 projects in 2020, and more than 3000 in 2050 representing an amount of 10 Gt of CO 2 stored per year. These figures highlight the large number of storage sites that shall be operated. Nevertheless, a prerequisite to any CCS large scale industrial development is the demonstration that the storage is both efficient and safe as outlined by the European directive on geological storage of carbon dioxide (directive 2009/31/EC, April 2009). In this context, monitoring techniques constitute a key component of any CO 2 storage risk management strategy aiming at providing information on the CO 2 fate (migration and trapping) within the reservoir and demonstrating that CO 2 containment is effective (e.g., [START_REF] Benson | Monitoring Carbon Dioxide Sequestration in Deep Geological Formations for Inventory Verification and Carbon Credits[END_REF].

Among the possible monitoring techniques, the space-borne Differential SAR (Synthetic Aperture Radar) interferometry may constitute a promising and cost-effective option, as it can provide valuable information on surface deformation induced by volumetric changes through time and space resulting from CO 2 injection in deep aquifers (Tamburini et al., 2010). This technique uses the phase of the radar signal in order to derive an estimation of the displacement of the ground surface between two passes of the satellite over the area. Further details can be found for instance in [START_REF] Raucoules | Use of SAR interferometry for detecting and assessing ground subsidence[END_REF]. The basic idea is that (once removed flat earth and topographic components) the phase difference is proportional (modulo 2) to the Line of Sight displacement occurred between the two acquisitions (assuming that we can identify a stable reference area on the image. If not, the displacement is obtained with an unknown constant value offset). A major interest of the technique lies on the fact that the spatial sensor is able to image the ground surface regularly (presently, repetitiveness varies depending on characteristics of the available satellites between few days to more than one month).

Yet, one major limitation results from the temporal decorrelation phenomena [START_REF] Zebker | Decorrelation in interferometric radar echoes[END_REF] stemming from changes in the physical (e.g. dielectric characteristics) or geometric (e.g. the distribution inside of the image elementary cell) characteristics of the backscattering objects. Multi-temporal processing of interferometric series is generally needed to cope with this shortcoming. Among the possible multi-temporal techniques, PSI (Persistent Scatterers Interferometrye.g. Ferreti et al., 2000) are advanced techniques based on the processing of numerous, typically more than 30, images to identify a network of persistent, temporally stable, highly reflective ground features -persistent scatterers (denoted PS). These scatterers typically are cultural features of the developed landscape such as buildings, utility poles, roadways, etc., but also natural stable elements with stable and strong backscattering characteristics (such as stable rocks with specific orientation).

To date, such space-borne techniques have been successfully used in a variety of domains involving subsurface fluid withdrawal / injection and surface deformation response, for instance groundwater flow characterization in complex aquifers systems [START_REF] Bell | Permanent scatterer InSAR reveals seasonal and long-term aquifer-system response to groundwater pumping and artificial recharge[END_REF], oil / gas fields characterization (e.g., Vasco et al., 2008a), monitoring of seasonal gas storage (e.g., [START_REF] Teatini | Geomechanical response to seasonal gas storage in depleted reservoirs: A case study in the Po River basin, Italy[END_REF]. More recently, surface deformation have been measured at In Salah Gas Project in Algeria, where 0.5-1 million tonnes CO 2 per year have been injected over the past 5 years (e.g. see site context described by [START_REF] Rutqvist | Coupled reservoir-geomechanical analysis of CO2 injection and ground deformations at In Salah, Algeria[END_REF]. This analysis revealed an average ground uplift of 5 mm/year and allowed getting better insight not only in the CO 2 migration at depth (e.g., [START_REF] Mathieson | Satellite imaging to monitor CO2 movement at Krechba, Algeria[END_REF], but also in the behaviour of the storage reservoir (e.g., Vasco et al., 2008b) and of the storage complex during injection operations (e.g., [START_REF] Rutqvist | Coupled reservoir-geomechanical analysis of CO2 injection and ground deformations at In Salah, Algeria[END_REF].

However, such a successful result might not be easily generalised to any onshore CO 2 storage sites, because the In Salah case is characterized by very favourable conditions for the implementation of the PSI technique, the site being located in an arid region so that the Permanent Scatterer PS density is very high, exceeding 1,000 per km² (300,000 PS used by Tamburini et al., 2010 over an area of approximately 10 km by 15 km).

The objective of the present paper is to get a better insight in the applicability of the Persistent Scatterers Interferometry implementation regarding the conditions of the future storage sites and the specificities of the expected ground deformation phenomenon, i.e. the expected spatio-temporal evolution of surface displacements caused by CO 2 injection operations at depth.

The remainder of the paper is organized as follows. We estimate in a first section the expected surface displacements by means of an analytical model. Though inherently based on simplifications, the analytical results allow investigating a large variety of injection scenarios (e.g. depth, injection rate, storage aquifer characteristics, etc.) useful to support the discussion in section 2. Given the discussed limitations for the CO 2 storage context, we propose to explore, in a third section, a possible strategy to enhance the performance of the PSI technique for CO 2 storage monitoring by deploying a network of "artificial" PS, i.e. corner reflectors CR (e.g., [START_REF] Novali | PSInSAR Validation by Means of a Blind Experiment Using Dihedral Reflectors[END_REF] in an iterative design approach over the whole injection duration.

Estimates of the expected spatio-temporal evolution ground vertical displacements

In this section, we briefly describe an analytical solution for assessing ground displacement assessment during CO 2 injection in deep aquifers (section 1.1). Details of the analytical model are provided in Appendix A and B. This model is then applied to a typical injection scenario (section 1.2) and its applicability is discussed (section 1.3).

Analytical solution for ground displacement assessment

In the present article, the expected spatio-temporal evolution of ground vertical displacements due to CO 2 injection at depth are estimated using an analytical model developed by combining :

-the approximate solutions for pressure build-up (i.e. the pressure increase due to injection) during CO 2 injection in brine aquifers of very large spatial extent (seen as open systems) based on the studies of Mathias et al. (2009a&b); -the poro-elastic solution for vertical surface displacements estimate in an elastic halfspace, originally developed in the field of hydrocarbon and geothermal activities by considering an axi-symmetric thin disk reservoir with general reservoir pressure distribution [START_REF] Segall | Induced Stresses due to Fluid Extraction from Axisymmetric Reservoirs[END_REF][START_REF] Segall | A note on induced stress changes in hydrocarbon and geothermal reservoirs[END_REF].

For sake of clarity, the main assumptions and formulations of both analytical solutions are respectively described in Appendix A and B. Noteworthy, each of these solutions have been validated through numerical simulations (Mathias et al., 2009a) or against in situ observations (at Lacq in southern France, see [START_REF] Segall | Induced Stresses due to Fluid Extraction from Axisymmetric Reservoirs[END_REF].

Analysis of a typical injection scenario

We assess the spatio-temporal evolution of ground vertical displacements resulting from a typical scenario corresponding to 1 Mt/y CO 2 injection during 25 years into a saline aquifer formation at 1500m depth. Considering a pressure and geothermal gradient of respectively 1.10 5 Pa/m and 0.03 °C/m, the initial pressure and temperature are of 15 MPa and 45 °C at 1500m depth. The homogeneous properties of the aquifer formation are typical of sedimentary formations suitable for CO 2 storage, with high-enough permeability and porosity (see for instance [START_REF] Birkholzer | Large-scale impact of CO 2 storage in deep saline aquifers: A sensitivity study on pressure response in stratified systems[END_REF] for typical hydrogeological properties and Rohmer and Seyedi, 2010 for typical geomechanical properties) and are summarized in Table 1.

[Table 1 about here] Figure 1 depicts the spatial distribution of changes in vertical ground displacements estimated using the equations described in Appendix A and B. Different injection durations are considered and vary from 1 month to 25 years. The pressure build-up reaches a maximum value within the injector near-zone and decreases deeply in an exponential fashion within the CO 2 plume region (the extent being marked by the steep decrease of the pressure build-up in Figure 1A) and linearly within the brine-filled region. The vertical surface displacements follow the same trend (Figure 1B), but with a smoother form of the spatial distribution so that these are approximately constant over approximately 1 km and steeply decrease reaching (and even exceeding) the lateral extent of the pressurized zone at depth (i.e. the zone for which the pressure build-up is non-nil). For instance, after 2 years of injection, the vertical surface displacements exceed 4 mm over the lateral distance range of 0 to 1 km and decrease reaching ~10 km given that the extent of the pressurized zone reaches ~8 km. The smooth spatial distribution is explained by cumulative contribution of all the rings of dilation with magnitude proportional to the pressure build-up (i.e. decreasing with the lateral distance) as shown in Appendix (Figure 9).

Let us consider the gradient (in absolute value) of displacements corresponding to the difference between vertical displacements respectively observed at positions r  and of

Δ    r divided by the distance Δ 
, where Δ  is the characteristic length of a satellite image pixel (typically between 10 and 50 m depending on sensor and processing). This gradient is largely smaller than the maximum allowed value corresponding to half a fringe per pixel, i.e. of about 1.5 cm/10m=1.5x10 -3 to 1.5cm/50m=3.10 -4 for vertical movements with ERS data [START_REF] Raucoules | Use of SAR interferometry for detecting and assessing ground subsidence[END_REF]. This value would be modified for sensors with different resolutions and wavelengths, the maximum measurable gradient being proportional to wavelength and inversely proportional to pixel size. Even in case of a point-based method (such as PSI techniques) with larger Δ   typically ranging 500m to 1000m in rural zones (in this case Δ  does not correspond to the pixel size, but to the distance between PS), this maximum measurable gradient should not be theoretically a limitation, because the expected total deformation of the order of ~1 cm, as shown in Figure 1B, corresponds to less than one Cband sensor fringe (equivalent to ~3 cm of vertical displacement).

[Figure 1 about here]

Considering the temporal evolution of the maximal surface displacement within the injector near-zone, Figure 2 shows that it reaches a maximum value of ~11 mm over a time span of 25 years and evolves following a logarithmical function over time. This trend is partly explained by the temporal evolution of the maximum pressure build-up within the injector near-zone.

But, given its rapid evolution (after 1 month of injection, it reaches more than 85 % of the maximal pressure build-up), the temporal evolution is in great part explained by the cumulative contributions of rings of dilation, as shown in Figure 9, (Appendix A), whose number increases over time as the pressure front propagates within the aquifer formation.

[Figure 2 about here]

Applicability of the analytical solution

The described analytical model is proposed as a quick-assessment tool, which should be used to capture the general trend of surface deformation due to CO 2 injection, i.e. to estimate orders of magnitude of the expected vertical displacements either to assess the likelihood of a risk event associated with ground surface perturbations in a preliminary risk analysis of the storage site (e.g., [START_REF] Bouc | Safety criteria for CO 2 geological storage: Determination workflow and application in the context of the Paris Basin[END_REF] measured.

An in-depth understanding of the surface deformation process due to CO 2 injection at depth should consider several aspects, which the analytical model hardly accounts for, most of them being related to the complexity of both the deformation processes and the geological subsurface structure.

The described model assumes poro-elastic rock behaviour, hence implying that uplift instantaneously follows linearly the reservoir volumetric increase. Though confirmed through in situ compaction measurements (see e.g. the Groningen field, studied by [START_REF] Mobach | In situ reservoir compaction monitoring in the Groningen field[END_REF], a nonlinear dependence may also appear in some cases. See detailed discussion by [START_REF] Hettema | Subsidence delay: Field observations and analysis[END_REF]. Furthermore, the ground deformation processes may be complicated by the presence of soils and unconsolidated rocks near the surface (e.g., [START_REF] Singh | Mine subsidence. SME Mining Engineering Handbook[END_REF].

The proposed model considers an idealized view of the subsurface structure so that both the reservoir and the surrounding rock are considered homogeneous and having the same material

properties. Yet, aquifer formations targeted for CO 2 storage are mostly located within large sedimentary basins, which can be viewed as multi-layered systems (see e.g., [START_REF] Birkholzer | Large-scale impact of CO 2 storage in deep saline aquifers: A sensitivity study on pressure response in stratified systems[END_REF]. Analytical models have been proposed to account for this aspect [START_REF] Fokker | Semi-analytic modeling of subsidence[END_REF] and can be extended as a future direction of research to the CO 2 storage case in the similar approach as proposed in section 1. Furthermore, the presence of geologic discontinuities (i.e. the existence of faults) may increase ground movement potential, but on a localized spatial scale, as numerically demonstrated by [START_REF] Rutqvist | Coupled reservoir-geomechanical analysis of CO2 injection and ground deformations at In Salah, Algeria[END_REF].

Constraints on the implementation of PSI technique

In turn, we examine constraints to the implementation of the PSI technique regarding the specificities of the deformation phenomenon associated with the CO 2 injection at depth, in terms of magnitude (section 2.1), spatial and temporal scales (section 2.2 and 2.3) and evolution (section 2.4). Table 2 summarizes the constraints on the PSI implementation and the needs in terms of sensor characteristics in order to address the issues resulting from the expected motion characteristics.

[Table 2 about here]

Magnitude of the deformation

It is difficult to determine beforehand whether the specific deformation phenomenon associated with the CO 2 injection at depth (as illustrated by the typical injection scenario described in section 1.3), may be monitored, because the feasibility and accuracy of PSI strongly depend on the number of satellite acquisitions and the temporal baselines (e.g. [START_REF] Hanssen | Satellite radar interferometry for deformation monitoring: a priori assessment of feasibility and accuracy[END_REF].

The value of the precision PSI techniques has been widely discussed in the literature and can vary from less than 1 mm to centimeters. [START_REF] Colesanti | Permanent Scatterers: Precision Assessment and Multi-platform Analysis[END_REF] suggested that in an operational context under very good conditions (large data setsi.e. >30 images, dense PS network, quasi linear deformation regimes), the precision could reach 0.5 mm/y for the annual deformation rate and 3 mm for the displacement values at given acquisition dates for PSI based on data from ERS, Envisat/ASAR or equivalent (C-band, resolution ~10m, repetitiveness of ~ 1 month). However, in less favourable cases (deformation rates higher than some cm/y, sparse PS, non-linear deformation) results could be worse. Studies based on intercomparison of different PSI techniques and validation respect to ground measurements (see for instance [START_REF] Raucoules | Validation and Intercomparison of Persistent Scatterers Interferometry: PSIC4 project results[END_REF] showed that in unfavorable cases, the precision value can drop to 1.5 mm/y and the accuracy can be much worse (>5 mm/y). Furthermore, deformation superior to ~cm/y tends to be under-estimated if no prior assumption on the deformation is used in the processing [START_REF] Raucoules | Validation and Intercomparison of Persistent Scatterers Interferometry: PSIC4 project results[END_REF]. This issue is not expected for CO 2 storage applications as shown by the estimates of ground displacements resulting from a typical CO 2 injection scenario (see Figure 1B).

On this basis, a detection threshold for an individual deformation value of 5 mm (i.e. below which the deformation phenomenon is very unlikely to be detected) is considered in the present study. For the considered injection case in section 1.3, we show that deformation phenomenon associated with the CO 2 injection is likely to be detected only after the first year of injection. Different injection scenarios are then considered with injection rates and depths respectively ranging from 0.5 to 2.5 Mt/y and 1000 to 2500 m and assuming the aquifer formation properties constant as described in Table 1. Figure 3A) shows the maximum vertical displacement reached at surface in the injector near-zone after respectively 1 and 10 years of injection (top and bottom left panels). Several aspects can be outlined:

 The maximal vertical displacements is proportional to the mass injection rates, whereas it is non-linear with depth mainly due to the non-linear form of the relationship between the fluid properties of CO 2 and brine with depth;

 The depth of the aquifer storage formation by 500m implies increasing the necessary injection rate by 0.5 Mt/y for the surface deformation to be detectable. For instance, ground displacements after 1 year caused by CO 2 injection at 0.5Mt/y within a 1000m

deep aquifer is equivalently detectable as those caused by CO 2 injection at 1.0 Mt/y within a 1500m deep aquifer;

 For the whole range of depths and time durations, the deformation associated with the 0.5 Mt/y injection scenario is hardly detectable ;

 Only the deformation associated with the 1000m deep injection may be detectable over the whole range of injection rates.

[Figure 3 about here]

Spatial scale

The basic idea of PSI techniques is that there exist points with stable backscatter characteristics during the period covered by the data series. In other terms, those points (or PS) have high (respect to other pixels of the images) backscattered amplitude that does not change significantly from one image to other. They are generally considered as objects with extents smaller than the pixels where they are located. Such PS are little affected by geometric and temporal decorrelation and allow (using statistical criteria) to compensate the atmospheric phase screens in order to improve the precision of the measurement. In urban area, buildings and other man-made structures often act as PS due to their corner reflector-like scattering behavior and high radar reflectivity so that we can obtain a very high density (hundreds per km²) of PS, even exceeding thousands per km² in arid regions (e.g., [START_REF] Bell | Permanent scatterer InSAR reveals seasonal and long-term aquifer-system response to groundwater pumping and artificial recharge[END_REF].

However, the total area at ground surface to be affected by the deformation phenomenon caused by CO 2 injection at depth is expected to be very large in relation with the large overpressurized zone at depth. Several authors have recently shown that the latter can reach of the order of several to more than a few tens of kilometers (e.g. [START_REF] Birkholzer | Large-scale impact of CO 2 storage in deep saline aquifers: A sensitivity study on pressure response in stratified systems[END_REF]. The detectable zone of deformation at surface will then be characterized by a lateral extent of the order of a few kilometers in the early stage of injection (ranging between 0.5 and less than 5 km after 1 year of injection, Figure 3B, top) and of more than ten of kilometers (ranging between ~5 km and less than 15 km after 10 years of injection, Figure 3B, bottom). Such a large size of the monitored area implies that the size of the area covered by the satellite acquisitions (i.e. swath) should ideally be larger than 50 km. Yet, the larger the swath, the lower the resolution of the sensor. For instance, TerraSAR-X (mode scanSAR) presents a very large swath of 100km, but with a resolution of ~16m, which may be insufficient (see discussion below in section 2.3 and 2.4). Conversely, very high resolution sensors (~1m, such as TerraSAR-X or Cosmo-Skymed) offer insufficient coverage (~10km) and their use implies acquiring and processing a much larger amount of image data, which may become cost and resource intensive.

Furthermore, over such large areas, we can reasonably expect a vegetated cover (fields, forests, etc.) on a significant part of the storage site. For illustration purpose, the PICOREF sector in the Paris basin (of 3,455 km², see location in Figure 6), where the potential implementation of CO 2 storage sites have been investigated (e.g. Grataloup et al., 2008), is characterized by a land use composed of less than 5 % of artificial / urban areas and more than 95 % of natural (mainly forests) and agricultural areas. In such natural terrains, the absence of bright man-made structures might lead to sparse PS network (with density ~10 PS every km² to less than 1/km²). This constitutes a major obstacle to the phase unwrapping stage, where the solutions are directly dependent on the PS density and makes reliable estimation of deformation using PS techniques a challenging task.

The use of L-band data may contribute to improve the processing given that the vegetation is better penetrated using signal with short wavelength (e.g. [START_REF] Raucoules | Use of SAR interferometry for detecting and assessing ground subsidence[END_REF].

Unfortunately there is currently, to the authors' knowledge, (and probably for the next 2 decades) no planned satellite mission fulfilling the requirements in terms of resolution and/or acquisition repetitiveness for the proposed monitoring application as discussed below.

Temporal scales

CO 2 injection operations are typically planned to last between 25 and 50 years. This "long term" monitoring implies considering several issues.

A first difficulty is the availability of a space-borne sensor, whose mission durations are typically of the order of a ten of years (e.g. see [START_REF] Hanssen | Satellite radar interferometry for deformation monitoring: a priori assessment of feasibility and accuracy[END_REF], Table 2). Thus, any long term monitoring of the injection operations should be designed potentially considering a minimum of 2 to 3 satellites, whose characteristics (orientation of the orbital trajectory, center frequency, C-, L-band, etc.) might change from a satellite to another and not be necessarily suited for interferometric analysis. For instance, the currently-used ENVISAT mission, which was deeply modified in October 2010 (change in the orbital characteristics, affecting the possibility of carrying out reliable PSI processing after this event), will come to an end in the coming years. A transition period with another space-borne sensor should be accounted for in the design of long term monitoring plan. For instance, the successor of ENVISAT, namely This can be considered an order of magnitude of the life expectancy of a PS.

Though requiring further investigations, the effect of the land cover nature (urban / rural) is expected to have little influence on the life expectancy of PS. It is true that on vegetated area the temporal decorrelation (on none-PS points) is globally higher than on urban areas, but pure point target PS are little affected by temporal decorrelation (as they are identified by the temporal stability of their signal) whatever on urban and on non-urban areas.

The use of high resolution data is expected to improve the life expectancy of PS. As shown by [START_REF] Crosetto | Persistent Scatterer Interferometry: Potential, Limits and Initial Cand X-band Comparison[END_REF], higher PS densities are expected with higher resolution sensors.

Assuming an equivalent PS decay rate for different resolutions, the total number of PS is expected to be larger with high resolution sensors considering an isolated small area with high PS density located in a larger area depleted of PS (e.g. small groups of buildings in a larger agricultural zone). Following this reasoning, the probability of having at least a few PS lasting during the whole observation period is therefore higher with high resolution sensors (provided that a maximum of isolated building groups will keep some PS during the observation period). Yet, as mentioned in section 2.2, the high resolution requirement is limited by the large swath requirement. Current sensor TerraSAR-X (Strip Map Mode) and future Sentinel-1 mission should both fulfill such a trade-off, but with major differences in the data provision policy.

This latter issue should be carefully accounted for, because long term monitoring possibly implies constituting a large data base (depending on the acquisition and processing frequency;

which is further discussed in section 2.4). To remain cost-effective regarding ground-based geodetic measurement techniques (e.g., optical levelling, GPS, tiltmeters, etc.), a possible strategy can be similar to the provision policy associated to the Global Monitoring for Environment and Security (http://www.gmes.info/), which allows obtaining the data at low price as envisaged for the future Sentinel-1 mission. The interested reader can refer to [START_REF] Rucci | Sentinel 1 SAR interferometry applications: The outlook for sub millimeter measurements[END_REF] for additional discussions on the use of future Sentinel-1 and new X-band sensors for reservoir applications.

Spatio-temporal evolution

The monitoring of CO 2 storage sites has two major objectives (see e.g., [START_REF] Benson | Monitoring Carbon Dioxide Sequestration in Deep Geological Formations for Inventory Verification and Carbon Credits[END_REF]. A first objective, categorized as "verification", focuses on the monitoring of the expected behaviour of the storage complex, in particular detecting the location and quantifying the size of CO 2 plume within the reservoir and tracking the pressure front propagation. The second objective, categorized as "early-warning", focuses on the detection of any "abnormal behaviour", i.e.

deviations from the expected behaviour of the storage complex. This latter case corresponds to the occurrence of unwanted events such as the sudden (i.e. over a short time span) multimillimetric uplift localized on a fault zone, which may be related to the fault reactivation at depth (see e.g., [START_REF] Rutqvist | Coupled reservoir-geomechanical analysis of CO2 injection and ground deformations at In Salah, Algeria[END_REF].

Regarding the "verification" objective, the acquisition and processing rate should be chosen based on the simulations of the expected spatio-temporal evolution of the deformation phenomenon. In this view, section 1 provides useful orders of magnitude for a typical injection scenario. At a given location, Figure 2 shows that the deformation phenomenon is characterized by a highly non-linear evolution over time. This might limit the performance, because the PSI technique is generally better applied for rather conservative types of deformation (e.g. [START_REF] Hanssen | Satellite radar interferometry for deformation monitoring: a priori assessment of feasibility and accuracy[END_REF], e.g., characterized with (approximately) uniform deformation velocity relatively linear temporal evolution of deformation velocity or periodic (such as seasonal natural gas storage, [START_REF] Teatini | Geomechanical response to seasonal gas storage in depleted reservoirs: A case study in the Po River basin, Italy[END_REF]. To tackle this issue, an alternative would be to focus on the discrepancies of the measured displacement field respect to the proposed evolution model, i.e. deviation to the predictive model. By subtracting the modelbased simulation of the phase to the interferometric data, an important part of the (non-linear) temporal signal would be removed, facilitating the PSI processing (in particular by reducing the risk of unwrapping errors). Besides, the important information is the (unforeseen)

discrepancy respect to the model that is needed for adjusting the representation of the phenomenon. The proposed iterative procedure described in section 3.3 that revisits periodically the knowledge of the phenomenon, could partly deal with this issue.

Figure 1B provides information on the spatio-temporal evolution of the "deformation front"

given a typical injection scenario (see details in section 1.2). An individual measurement point located at a given lateral distance from the injection zone, requires a major and more frequently sampled monitoring during a specific period of the order of ~2-4 years, the beginning of this period depending on the distance to the injector. In other terms, a "deformation front", as schematically depicted in Figure 4, affects the location during this short time span t=2-4 years. Before the front has passed no deformation affects the point.

After, the point has smaller displacement (see Fig. 4). Such a spatio-temporal behaviour has implications on both the choice of the frequency of image data acquisition (i.e. sampling procedure) and of the processing rate (i.e. frequency of PSI analysis).

At a given location, the monitoring should be highly-rated sampled during the short period of the deformation front passage. As indicated in section 2.1, >30 images is an ideal data set size for precision reasons. The orbit cycle of ASAR-like sensors (about 10 acquisitions per year) may be therefore insufficient to carry out such a measurement at the limit of the sensitivity of the method. Current sensor TerraSAR-X (Strip Map Mode) and future Sentinel-1 mission should fulfill such a requirement (with both repetitiveness of the order of a ten of days), but with differences in the data provision policy (see section 2.3).

Besides, changes in the velocities values of the differential displacements are expected to be detectable only after a few months (except for the case of multi-millimetric sudden ground movements, which is discussed in the following). The processing of the data image with very high frequency (e.g. every day) may not be appropriate given the described physical evolution of the deformation, but also considering cost and resource constraints. Considering the expected spatio-temporal evolution of the deformation phenomenon as calculated for the typical injection scenario described in the section 1.2, one PSI analyse every 3 -6 months can be proposed to fulfil the "verification" objective of monitoring. Noteworthy, this processing rate is intermediate between quasi-real time (e.g. bottom-hole fluid pressure) and pluriannual monitoring techniques (e.g. seismic campaign).

Regarding the "earlywarning" objective of the monitoring, two cases can be considered. On the one hand, potential zones at risk (e.g. fault zones, injection zone, etc.) can be identified prior the injection operations based on a risk assessment procedure (see [START_REF] Bouc | Safety criteria for CO 2 geological storage: Determination workflow and application in the context of the Paris Basin[END_REF]. On these localized zones, a specific PSI analysis can be envisaged using high resolution/high repetitiveness (but reduced swath) data equivalent to those provided by the current X-band sensors TerraSAR-X or COSMO-SkyMed (see e.g., Ferretti et al. 2011b), to complement the PSI analysis for the "verification" monitoring. Noteworthy, other ground-based geodetic techniques would require systematic on-site measurements, such an approach being time consuming and resource intensive.

On the other hand, if any "abnormal behavior" is detected in an "un-expected" zone either after processing the image data or based on other indicators (induced seismicity recorded at the surface due to possible fault reactivation; increase of the fluid pressure recorded in the monitoring aquifer indicating a possible leakage, etc.), the considerations for data acquisition and processing rate as discussed for the verification monitoring should be adjusted. The great advantage of PSI analysis compared to ground-based geodetic measurement is the possibility to assess deformation after the unwanted event has occurred, because image data would have been archived without user's intervention (e.g., every12 days using the future Sentinel-1 sensor). Punctual and local studies on short time spans can also be proposed to focus on the "un-expected" zone at risk to complement the PSI analysis for the verification monitoring.

[Figure 4 about here]

Planning a network of artificial reflectors

Description of corner reflector devices

One major difficulty to the implementation of PSI in the CO 2 storage context stems from the potentially low PS density over such large areas as previously discussed. Though improved co-registration and use of more efficient computational techniques (the interested reader can refer to [START_REF] Sousa | Persistent Scatterer InSAR: A comparison of methodologies based on a model of temporal deformation vs. spatial correlation selection criteria[END_REF] can provide valuable support, a potential solution to enhance the performance in regions where bright man-made structures are merely absent can rely on the installation of artificial PS, i.e. corner reflectors CR (e.g., [START_REF] Novali | PSInSAR Validation by Means of a Blind Experiment Using Dihedral Reflectors[END_REF]. These artificial metallic devices (of generally trihedral or cubic geometry characterized by a few tens of centimeters long) are installed on ground providing high backscatter to the radar signal to complement the existing "natural" PS network, i.e. to fill the gap between natural persistent reflectors. This approach has been proposed in other contexts (mine, volcano, coastland and marshland environments, etc.), where the ground movements' detection has proved to be very challenging. See for instance the recent application case in the Venice lagoon [START_REF] Teatini | Resolving land subsidence within the Venice Lagoon by persistent scatterer SAR interferometry[END_REF].

Optimizing the position of the CR

To the authors' knowledge, very few works have handled the issue of planning the position of CR in a predictive fashion, i.e. before any ground deformation measurement has been carried out. For instance, in the recent interferometric analysis of the land subsidence in the Venice lagoon [START_REF] Teatini | Resolving land subsidence within the Venice Lagoon by persistent scatterer SAR interferometry[END_REF], a network of 55 CR has been planned taking into account a fixed distance of 1 km between them of any adjacent natural reflectors. Yet, these were installed only once the full interferometric analysis had been carried out for a ten of years.

Planning the CR network using an "one-size-fits-all" approach, i.e. using a fixed distance value may rapidly become impracticable considering the large area expected at the end of injection (of ~325 km² after 25 years for the case modeled in section 1.3). Optimizing the location of the CR should result from a trade-off between the issues related to the installation of such devices (costs, authorizations, maintenance, etc.) and the requirements in terms of coverage and precision for the PSI results. Furthermore, the time-scales of CO 2 storage constitute a major constraint for the monitoring plan that should be designed to follow the long term evolution of the ground movement over a time span ranging from 25 to 50 years.

An example of a ground-based network, i.e. of devices fixed at the surface that exist over such a long time duration, is provided by the networks of benchmarks used for national spatial reference systems (e.g. in France, the currently-used network has been installed since the 19 th century), but dimensions of such ground-based devices have nothing comparable to the ones of a CR.

In this context, a crude design approach based on the expected surface deformation calculated regarding the whole time duration of the storage project may be highly limited by both maintenance issues of the CR network given such large spatio-temporal scales and by the changes in the sensors characteristics (typically each 7-10 years, see section 2.2) over the full lifecycle of the storage injection phase. Therefore, we advocate a design strategy thought in an iterative fashion: [START_REF] Fokker | Semi-analytic modeling of subsidence[END_REF] or in a first approach, the one described in section 1 presenting the attractive feature of a rapid estimation, but limited by the simplifications inherent to the analytical nature of the model of such a simplified model (as discussed in section 1.3). A more sophisticated analysis can rely on large-scale coupled numerical models as demonstrated by [START_REF] Rutqvist | Coupled reservoir-geomechanical analysis of CO2 injection and ground deformations at In Salah, Algeria[END_REF] in the In Salah's context. The advantage of such analysis are the ability to take into account complex injection scenarios, complex reservoir pore pressure spatio-temporal evolution, complex geological structures (e.g., multi-layered system as present in sedimentary basin suitable for CO 2 storage, [START_REF] Birkholzer | Large-scale impact of CO 2 storage in deep saline aquifers: A sensitivity study on pressure response in stratified systems[END_REF]  Step2: based on the expected extent of detectable deformation surface, localize the a priori good candidate for PS, i.e. spatial points that appear to be persistent in the present, but above all, that shall remain persistent over the full duration of the injection. Long term "natural" PS can constitute man-made structures within urban zones (historical monuments, churches, industrial plants, roads, etc.). We can consider that a large number of these spatial points identified as PS will keep their characteristics (reflectivity) over the time duration of the considered epoch (of a ten of years). This assumption can be considered valid regarding the study of Ferreti et al.

(2003), as discussed in section 2.2. ;

 Step3: optimize the position of the corner reflectors CR considering the expected detectable surface deformation for the whole duration of the injection. We advocate choosing the distance between CR and natural PS based on a criterion dependent on the initial network of natural PS. The idea is to maximize the spatial coverage of the area where ground deformation is expected to be detected considering the location of identified potential "natural" PS. Formally, this can be achieved by treating the "natural" persistent scatters as Voronoi sites and by affecting to each of them a Dirichlet cell (also named Voronoi cell) consisting of all points closer to PS than to any other site as schematically depicted in Figure 6A. The segments of Voronoi cell are all the points in the space that are equidistant to the two nearest sites PS. The potential positions of the CR are then chosen as the Voronoi nodes corresponding to the equidistant to three (or more) PS sites. The positions of the CR merely derive from a pure mathematical space-filling criterion, but in real-case applications, adjustments will certainly be necessary if the Voronoi-based method indicates installing the CR in densely vegetated cover (forests, wetlands, etc.), which may highly prevent the detection or the installation of the CR, or too close to another CR (e.g., with a distance <1 km). Then, the CR only lying within the expected detectable zone for the considered 10-years epoch, can be selected. These CR should be installed to help monitoring the expected surface deformation over the duration of one epoch.

Based on the discussions in section 2, the PSI analysis to fulfil the "verification" objective of the monitoring can conducted using data characterized by a large swath, a rather good resolution and repetitiveness of a ten of days. The processing rate of the image data should be chosen depending on the simulation predictions for the spatio-temporal evolution of the deformation phenomenon associated with CO 2 injection. The simulations for the typical injections scenario (section 1.2) show that a processing of 3-6 months can be envisaged during the 10-year epoch. This rate can be adjusted especially depending on the detection of a sudden "abnormal behaviour" on a localized zone (e.g., sudden multi-millimetric uplift due to possible fault reactivation). A complementary PSI analysis using high resolution / high repetitiveness (but reduced swath) data can be proposed on these zones at risk (see discussion in section 2.4). Such an iterative approach allows a frequent re-interpretation of the discrepancies between model and measurements in terms of evolution predictions (i.e.

deviation to the predictive model). Valuable information on the storage complex behaviour can be extracted in an inversion procedure as proposed more specifically in the field of CO 2 storage by Vasco et al. (2008a&b). The combination of both information (correction of the model based on the measurements and inversion of the complex behavior) can serve to increase the quality and robustness of the predictions.

The procedure is then repeated for the next 10-years epochs (step 1). Using the initial PS network enriched with the CR identified for the previous epoch (step 2), the next CR (i.e.

necessary for monitoring the next epoch) can be identified using the updated expected detectable surface deformation region (step 3), and the procedure is then repeated.

Noteworthy, Before the end of each 10-years observation epoch, a transition period should be envisaged between two successive satellite missions, with sensors of possibly different characteristics. During this period, a sufficient amount of SAR data (> 30 images) should be acquired for the PSI baseline using the new sensor. Considering a repetitiveness of a ten of days, this period should last ~1 year. Figure 5 summarizes the sequence of the steps considering two epochs and two different sensors.

[Figure 5 about here]

In case of sensors with different characteristics, a change in the configurations (e.g. orientation) of the CR, which have been installed during the first epoch, should be envisaged to make these CR "visible" by the new sensor. In such a case, a gap in the produced deformation time series may be expected of the order of few months, but only partly on the ground surface covered by the SAR data, i.e. corresponding to the region whose CR were adapted to the first 10-years epoch. To the authors' knowledge, the transition between two sensors during the implementation of SAR interferometry combined with CR has not been experienced yet and requires further investigations.

Finally, considering such a long term monitoring, the possibility that during certain periods, the lack of satellite mission with suitable characteristics cannot be fully discarded. In such a case, the information provided before and after the information gap could still provide useful information on the evolution of the injection processes.

Application to a potential storage site in the Paris basin

To illustrate the main steps of the proposed methodology, we apply the described methodology to a potential injection site selected within the PICOREF sector in the Paris basin (Figure 6). This potential storage site was identified based on the combination of criteria in terms of capacities and injectivity, environmental constraints, existing land-use and underground-use, economic and social aspects and based on a comprehensive risks analysis (Grataloup et al., 2008;[START_REF] Bouc | Safety criteria for CO 2 geological storage: Determination workflow and application in the context of the Paris Basin[END_REF]. Figure 6 shows the location of the potential injection well (outlined by a red cross marker). We consider 40 years of CO 2 injection at a rate of 2 Mt/y into the 1500m deep carbonate aquifer of the Dogger formation. The aquifer properties are: intrinsic permeability of 1.10 -13 m²; porosity of 12 %; Young's modulus of 24

GPa; Poisson's ratio of 0.29 (Rohmer and Seyedi, 2010 and references therein).

[Figure 6 about here]



Step 1: in this example, the analytical model described in section 1 is used to estimate the lateral extent of the expected detectable zone of vertical displacements  Step 2: we consider each centroids of the urban zone present within the detectable zone for the whole injection duration (total number of 24 zones). These represent localized zones where good candidates for natural PS are highly likely to be detected.

In a preliminary approach, we only consider these centroids as natural PS (outlined in Figure 7 as pink diamond markers);

 Step 3: on this basis, we assign to each of the centroids a Voronoi cell (Figure 7A) and choose the potential CR as the Voronoi cells (outlined in Figure 7A as blue triangular markers); the CR lying within forest zones or distant from each other of less than 1km are detected and are re-adjusted (deleted or re-localized). A total number of 23 CR are in this manner identified within the detectable region for the whole duration of the injection (here 40 years, black -coloured circle in Figure 7B).

[Figure 7 about here]

Using the region of detectable deformation for the first 10 years of injection (light greycoloured circle in Figure 7B), a total of 6 CR are identified through the Voronoi analysis.

Among the identified CR, two CR should be re-localized, because they lie within forest areas and one CR should deleted, because it is too close from the others. Without any other consideration (such as legal issues, accessibility of the site, …), we are able to propose a monitoring procedure adapted to the expected characteristics of the deformation for the first 10 years of injection both based on CR-InSAR and on current PSI techniques.

For the next 10-years epoch, Figure 8 provides an illustration of the location of the corresponding CR (i.e. after 10 years of injection) using the original PS network (Figure 7A, pink coloured diamond markers) enriched with the 5 CR selected for the first 10-years epoch.

The expected detectable regions for the whole injection duration (at 20, 30 and 40 years) are then re-estimated (step 1), here assumed unchanged for illustration purposes. On this basis, the CR for the whole duration are identified (step 2 & 3). Then, using the expected detectable deformation region for the next 10-years epoch (here for the date 20 years), the CR for this epoch are selected (Figure 8B).

[Figure 8 about here]

Concluding remarks and recommendations

The present paper proposes a description of the characteristics of the surface deformation, which will be expected during the injection of CO 2 into deep aquifers during 25 -50 years.

On this basis, the associated constraints for the planning of a monitoring plan relying on PSI techniques are discussed. To overcome such limitations, we focused on the optimisation of a CR (corner reflectors) network. An iterative strategy is proposed, which is justified by the large spatial extent of the expected detectable region and by the changes in the sensors characteristics (typically each 7-10 years) that could result in a need of re-organisation of the network (and even CR characteristics such as orientation). At each period (say, every 10 years), the number and location of the CR, necessary for the monitoring of the surface deformation of the next epoch, are estimated based on the prediction of surface deformation using large-scale geomechanical models. Furthermore, we propose to focus on discrepancies between model and measurements, which can be useful in terms of evolution predictions. To fulfil the "verification" objective of the monitoring, the PSI analysis should be frequently conducted during each epoch with a processing rate, which should be chosen depending on the predicted spatio-temporal evolution of the deformation phenomenon. The simulation results for a typical injection scenario show that a processing rate of the order of a few months may be sufficient. Such a frequent re-examination of the situation and model correction is therefore useful in a perspective of forecasting consequences of the injection. To fulfil the "early-warning" objective of the monitoring, a complement should be envisaged so that well definite zones at risk (such as fault zones) should be monitored using high repetitiveness / high resolution (but reduced swath) data. In both cases, a careful attention should be paid to cost and resource issues especially regarding provision policies of image data.

Finally, from now to the start of injection and during the 50 years injection, improvement of PSI processing techniques can occur. For instance, Ferretti et al. (2011a), recently proposed a technique (SqueeSAR TM ) based on the concept of diffuse scatterers, which can improve considerably the density of usable "natural" scatterers on certain un-urbanised area. See the recent application of such a technique at In-Salah [START_REF] Rucci | Sentinel 1 SAR interferometry applications: The outlook for sub millimeter measurements[END_REF]. If such improvement is efficient for surface characteristics such as those existing on the monitored area, an optimisation of the network could be proposedprobably based on a reduced CR number.

The periodic re-thought of the CR network will therefore allow the possibility of methodological improvement in addition to taking into account new sensors characteristics. The acquisition and processing rate should be adapted to the expected spatio-temporal evolution of the deformation phenomenon. Based on simulations on a typical injection scenario:

Table caption

 the regions, where the deformation front passes, should be highly-rated sampled for a 2-4 years period;

 the changes in the velocities values of the differential displacements are expected to be detectable only after a few months, implying a processing rate of 3 to 6 months.

Potentially sudden unexpected ground movements on localized areas

To fulfil the "early-warning" objective of the monitoring, a complement should be envisaged so that well definite zones at risk (such as fault zones) should be monitored using high repetitiveness / high resolution (but reduced swath) image data.

Figure captions

  Sentinel 1, is to be launched in 2013. This 7-12 years mission will consist in a space-borne SAR in C-band (and therefore providing data similar to ASAR) with a repetitiveness of 12 days. Considering an ideal data-set for PSI baseline of >30 images the transition period should last ~10 months, i.e. before the planned commercial-scale deployment of storage sites by 2020(IEA, 2009).A second difficulty arises from the decorrelation of PS over time, i.e. the "life expectancy" of the PS over such long time durations.Ferreti et al. (2003) tackled this issue by identifying points corresponding to PS only on sub-sets of the full ESA-ERS data set over the period 1992-2003 on the Milano urban area. This consisted in identifying PS, which have appeared or have disappeared during a given time span between the first and last acquisition. The number of such points is usually 10% of the number of PS over a period of a ten of years.



  Step1: estimate the expected spatio-temporal evolution of ground surface displacements induced by CO 2 injection operations at different epochs for the full period of injection. The duration of each epoch should be no longer than a ten of years to handle the possible change in the sensors characteristics. Let us consider the estimates of the surface deformation after 5 epochs of 10 years each. This step should be based on large-scale geomechanical analysis of the complex site. This can rely on analytical models such as the multilayered visco-elastic solution developed by Fokker

(

  considering a 5 mm detection threshold) induced by a whole CO 2 injection duration of 40 years and considering four 10-years long epochs. The final zone represented by a black circle in Figure 6 covers a total area of ~220 km² and corresponds to the zone of maximum extent that should be monitored by the PSI technique over the whole life span of the potential storage site;

Figure 1 :

 1 Figure 1: A) Pressure build-up (Pa) within the storage aquifer for a 1Mt/y CO 2 injection

Figure 2 :

 2 Figure 2: Temporal evolution of the maximum vertical displacement u max (mm) within the

Figure 3 :

 3 Figure 3: A) Maximum vertical displacement reached at the surface respectively after 1 and

Figure 4 :

 4 Figure 4: Schematic representation of the deformation to be measured. t corresponds to the

Figure 5 :

 5 Figure5: main steps for the application of the PSI analysis for "verification" monitoring using

Figure 6 :

 6 Figure 6: Location of the PICOREF sector (Grataloup et al. 2008) in the Paris basin context

Figure 7 :

 7 Figure 7: A) Voronoi analysis: pink coloured diamond markers represent natural PS located

Figure 8 :

 8 Figure 8: A) network of natural PS (pink coloured diamond markers) completed with the CR

Figure 9 :

 9 Figure 9: Schematic representation of the cumulative contribution (light blue coloured arrows)

Table 1 :

 1 Aquifer homogeneous properties used for the spatio-temporal analysis of the surface

	vertical displacements
	Aquifer property	Symbol	Value	Unit
	Intrinsic permeability	k	1.e -13	m²
	Porosity		15	%
	Thickness	H	50	m
	Salinity	sal	150	g/l
	Rock matrix	c rock	4.5e -10	Pa -1
	compressibility
	Young's modulus	E	15	GPa
	Poisson's ratio		0.25	-
	Biot's coefficient	b	1.0	-

Table 2 :

 2 Summary of the constraints on the implementation of the PSI technique and on the sensor regarding the specificities of the deformation phenomenon associated with CO 2

	injection at depth
	Expected deformation	Corresponding constraints on the implementation of the PSI
	characteristics	technique and on the sensor
	Large spatial scales tens	 Large swaths (>50 km) are required;
	of km)
	Large vegetated areas	 Poor PS density is expected. Additional Corner Reflectors
		could be required to complement "natural" PS network.
	Long term monitoring
	(25 -50 years)

 Because of limited sensors lifetimes, the use of different successive satellites is required;  High sensor repetitiveness (typically, ten of days or better)

Acknowledgements

This work was funded under the BRGM's Directorate of Research project. We also thank Marcello de Michele for useful discussions.

Appendix A: General solution to ground displacements Segall and co-authors [START_REF] Segall | Induced Stresses due to Fluid Extraction from Axisymmetric Reservoirs[END_REF][START_REF] Segall | A note on induced stress changes in hydrocarbon and geothermal reservoirs[END_REF] have presented poro-elastic solutions for stress and displacement fields estimate in an elastic half-space, induced by fluid extraction (or injection, the solutions being merely the opposite) for general axi-symmetric geometries with arbitrary pressure radial distribution within the reservoir (given as a function of space and time for instance from field data and/or reservoir simulations). Noteworthy this solution provides both horizontal and vertical surface displacements, but we will focus on the latter one.

The variation in reservoir pore pressure p is assumed to mainly govern the volumetric strain in the reservoir. The vertical displacement at the radial distance r at the free surface u z (r) results from the cumulative contribution of all the rings of dilatation at radius r ~ and depth Z (as schematically depicted in Figure 9), the magnitude of dilation being proportional to the product of the pressure p and of the Biot's coefficient b. Based on the Green's function at the free surface g z (r) for the radially symmetric case (Equation A1b), the solution for u z (r) holds as follows:

(A1)

where:

(A2)

where:

and F(K) is the complete elliptic integral of the second kind.

The physical properties are described in Table 1.

[Figure 9 about here]

Appendix B: General solution to the pressure build-up

The radial pore pressure p(r, t) evolution (Equation B1) is assessed based on the analytical solution developed by Mathias et al. (2009a) using the method of matched asymptotic expansions for the non-inertial flow of CO 2 into a slightly compressible brine aquifer (assumed to be an open system, and the pressure corresponds to the initial one) with a vertical pressure equilibrium assumption. This solution holds as follows, distinguishing three zones:

-the CO 2 fully saturated zone (near the injector), for x2 (x being dependent on r, see Eq. B2 and  is the ratio of the viscosity of CO 2 µ CO2 to the one of the brine µ brine );

-the diphasic zone (brine + CO 2 ), for 2  x  2/ ; -the brine-filled zone in the far field for 2/  x.

(B1)

where:

(B2)

where:

κ ≈ 0.5772 is Euler-Mascheroni constant ; r W is the injector well radius (m);  CO2 is the CO 2 density (kg/m 3 ); Q W is the mass injection rate (kg/s); t time (s); c brine is the brine compressibility (Pa -1 ); Other parameters are described in Table 1. Fluid properties are derived (viscosity, density and compressibility for both brine and CO 2 ) based on the empirical correlations with pressure and temperature supplied by Mathias et al. (2009b).


Franc ¸ois Keith

Nicolas Mansard
email: nmansard@laas.fr

Sylvain Miossec
email: sylvain.miossec@bourges.univ-orleans.fr

Abderrahmane Kheddar
email: kheddar@lirmm.fr

Optimization of Tasks Warping and Scheduling for Smooth Sequencing of Robotic Actions

This paper presents a method for sequencing a set of robotic tasks in an optimal way. Tasks description and execution are based on the task-function approach, which enables to build complex whole-body behaviors from local control laws. A naive solution to this problem would be to schedule the execution of the tasks sequentially, avoiding concurrency. This solution does not exploit full robot capabilities such as redundancy and have poor performance in terms of execution time or energy. However, reasoning on concurrent tasks is difficult while accounting for all the physical constraints of the robot. Our contribution is to determine the time-optimal realization of the mission taking into account robotic constraints that may be as complex as collision avoidance. Our approach achieves more than a simple scheduling; its originality lies in maintaining the task approach in the formulated optimization of the task sequencing problem. This theory is exemplified through a complete experiment on the real HRP-2 robot.

I. INTRODUCTION

A robot is designed to perform missions in various application contexts. When the environment is well or partially structured most missions can be hierarchically decomposed. That is, missions undergo functional objective decomposition into a set of processes or operations that can be defined as templates. Each operation can be decomposed into a set of tasks (i.e. generic sensory-motor functions), and each task can be easily mapped into robot execution. The whole scheme may constitute an exploitable generic skill/behavior. Yet, various levels of decomposition can be achieved depending on the envisaged software/hardware implementation, additional environment constraints, the human-machine interface, etc. In the end, the robot is assigned with a sequence of tasks to realize a given mission.

Numerous works have been proposed to compute such a sequence of tasks from a given mission and a set of causal paradigms [START_REF] Dechter | Constraint Processing, chapter 12, Temporal Constraint Network[END_REF], [START_REF] Ghallab | Automated Planning: Theory and Practice[END_REF]. However, they generally produce a symbolic plan, where the only numerical precisions lie on the scheduled time data. Its robotic application into the real world requires the time sequence to be refined, typically through an applicative path planner [START_REF] Lavalle | Planning Algorithms[END_REF], that will compute the trajectories to be followed by the robot. Yet, the meaning of the symbolic plan is lost in the global trajectory. Such low-level methods lack of robustness to environment changes or uncertainties. Consequently, the remaining trajectory may have to be recomputed several times while the mission is being achieved. Moreover, it is difficult (and then often specifically hard coded) to enhance the numerical trajectory with symbolic data, that would help re-computing only part of the plan [START_REF] Py | Dependable execution control for autonomous robots[END_REF] or distort locally the trajectory to apprehend small environment changes [START_REF] Quinlan | Elastic bands: Connecting path planning and robot control[END_REF], [START_REF] Lamiraux | Reactive path deformation for nonholonomic mobile robots[END_REF].

Rather than using a trajectory planner between the temporal reasoning and its real robotic execution, we propose to use a sensory-motor control approach based on task components. The task function [START_REF] Samson | Robot Control: the Task Function Approach[END_REF] or the operational space formulations [START_REF] Khatib | A unified approach for motion and force control of robot manipulators: The operational space formulation[END_REF] are elegant approaches to produce intuitively robot objectives. They also allow to address the control problem directly in the sensor space, improving robustness of the action execution against environment uncertainties and variability [START_REF] Espiau | A new approach to visual servoing in robotics[END_REF]. They are trajectory free, which means that it is not necessary to explicitly compute all trajectories before the execution or during the execution, namely in response to environment changes. Moreover, since a same task space is valid for a large set of robots, a control scheme based on task formalism is certainly portable and easy to modify and to maintain. In addition, these methods include a kinematics or dynamic formulation that decouples the task space from the null-space (i.e. the joint space that let the task invariant) [START_REF] Liégeois | Automatic supervisory control of the configuration and behavior of multibody mechanisms[END_REF], [START_REF] Hanafusa | Analysis and control of articulated robot with redundancy[END_REF]. A secondary task can then be applied in the null-space, and, recursively, a hierarchic set of tasks (or stack of tasks, or SoT) can be considered [START_REF] Nakamura | Task-priority based redundancy control of robot manipulators[END_REF], [START_REF] Siciliano | A general framework for managing multiple tasks in highly redundant robotic systems[END_REF]. Hierarchy of tasks are becoming popular to build complex behavior for very redundant robot such as humanoids [START_REF] Baerlocher | An inverse kinematic architecture enforcing an arbitrary number of strict priority levels[END_REF], [START_REF] Mansard | Visually-guided grasping while walking on a humanoid robot[END_REF], [START_REF] Sian | A switching command-based whole-body operation method for humanoid robots[END_REF], [START_REF] Sentis | A whole-body control framework for humanoids operating in human environments[END_REF]. The formalism introduced in [START_REF] Mansard | Visually-guided grasping while walking on a humanoid robot[END_REF] proved to be efficient in dealing with complex humanoid missions: the SoT is mainly a hierarchy of tasks driving the robot while ensuring locally a given set of constraints to be satisfied. We make use of this formalism (Section II).

A task (i.e. a task function [START_REF] Samson | Robot Control: the Task Function Approach[END_REF]) can be directly linked to the symbols on which the task temporal network is reasoning (e.g. reaching an object to be grasped is a task that requires the robot arm to be available, and whose post-condition is to have the gripper on the object -it is also directly described by a sensory-motor function applicable to the SoT).

Mission decomposition is thus executable directly by a SoT, which guarantees good robustness and avoid unnecessary trajectory (re)computation. However, exclusive task sequencing on the robot produces generally jerky movements which may look to humans as monotonous automated motions. On the other hand, it is difficult for the temporal network to produce a scheduling with task overlapping when the tasks concurrency is restricted by physical limitations of the robot (for example obstacles or dynamical constraints on a humanoid). Since the problem is not in a standard discreet form, symbol-based task scheduling techniques can not apply straightforwardly. On the other hand, semi-infinite optimization techniques [START_REF] Lee | Newton-type algorithms for dynamics-based robot movement optimization[END_REF] [START_REF] Miossec | Development of a software for motion optimization of robots-application to the kick motion of the HRP-2 robot[END_REF] have been used to generate low level trajectories for the overall execution, while accounting for such constraints. Such trajectory-based approaches are generally insufficiently robust to environment uncertainties.

In this paper, we propose to rely on task for both the overall symbolic reasoning and the control on the robot.

In between, we propose to use semi-infinite optimization to refine the symbolic schedule and account for (numeric) robotic constraints. Given a set of elementary tasks sequence to achieve a given mission, our solution returns for each task, the optimal times at which it is put and removed from the SoT and also the optimal parameters for the task execution (e.g. control gain). We additionally expect from this method a smooth tasks sequencing (i.e. smooth transitions of tasks through task overlapping). The originality of our approach lies in keeping the task component in the optimization formulation of this problem, which can roughly translate to optimizing tasks overlapping by manipulating tasks, i.e. the controllers, as variables of the optimization problem. The task formulation details are first recalled in Section II. A generic solution for optimizing a task sequence is then detailed in Section III. The theory is finally exemplified through an experiment with a real HRP-2 robot, the mission consisting in getting a can from a fridge.

II. GENERIC TASK SEQUENCING

A. Task function formalism and Stack of Tasks

Defining the motion of the robot in terms of task simply consists in choosing several control laws to be applied on a subpart of the robot degrees of freedom (DOF).

A task is defined by a vector e (typically, the error between a signal s and its desired value, e = ss * ). The Jacobian of the task is noted J = ∂e ∂q , where q is the robot configuration vector. In the following, we consider that the robot input control is the joint velocity q. The equation of motion is thus reduced to the kinematics:

ė = J q (1)
Considering a reference behavior ė * to be applied in the task space, the control law to be applied on the robot whole body is given by the least-square solution:

q = J + ė * + Pz (2)
where J + is the least-square inverse of J, P = I -J + J is the null-space of J and z is any secondary criterion that will be applied without disturbing the main task thanks to the projection into P 1 . A typical requested behavior is the 1 Eq. ( 2) is the least-square solution when z = 0 regulation of the error, which can be obtained through an exponential decrease by setting:

ė * = -λe (3)
As mentioned earlier, (2) enables to compose a complex behavior from a set of tasks [START_REF] Siciliano | A general framework for managing multiple tasks in highly redundant robotic systems[END_REF], [START_REF] Baerlocher | An inverse kinematic architecture enforcing an arbitrary number of strict priority levels[END_REF], [START_REF] Sentis | A whole-body control framework for humanoids operating in human environments[END_REF]: z can be used to fulfil a secondary task, without disturbing the main task having priority. This nice decoupling can be extended recursively to a set of n tasks, each new task being fulfilled without disturbing the previous ones. The complete implementation of this approach is proposed in [START_REF] Mansard | Task sequencing for sensor-based control[END_REF] under the name Stack of Tasks (SoT). The structure enables to easily add or remove a task, or to swap the priority order between two tasks, during the control. Constraints (such as joints limit) can be taken into account. The continuity of the control law is preserved even at the instant of change.

B. Gain handling

The simple attractor presented in (3) has the advantage to introduce a nice exponential decrease. However, it can be penalizing, since q is directly proportional to e (3) . At the beginning of the task, e reaches its higher value (strong acceleration), while at the end of the task, e decreases slowly (slow convergence).

A very classical 'trick' when regulating a task is to rather use an adaptive gain λ = λ(e(t)) that depends on the norm of the error of the task. To keep the nice property of the attraction, the gain only adapts with the error, and not directly with the time. We choose the following function:

λ(e) = (λ F -λ I ) exp -e β λ F -λ I + λ I (4)
with λ I the gain at infinity, λ F the gain at regulation (such as λ I ≤ λ F ) and β the slope at regulation.

C. Sequence of tasks

A task sequence is a finite set of tasks sorted by order of realization, and eventually linked to each other. Any pair of tasks can be either independent (i.e. they can be achieved in parallel) or constrained (i.e. one may have to wait for another one to be achieved, so as to make sense or to be doable).

The sequence can be formulated into a classical temporal network scheduling, starting at t 0 and ending at t End . Both values are finite and the sequence does not loop. Besides, we may consider for the sake of clarity but without loss of generality that each task appears only once in the sequence.

The position of a task in the sequence is defined by the time interval during which it is maintained in the SoT. For a given task i, this interval is noted [t I i , t F i ]: the task enters in the SoT at t I i and is removed at t F i . These instants are defined with respect to the beginning of the sequence at t 0 . However, they do not indicate the achievement level of the task: t F i may apply before the task regulation. Let's ǫ i be the tolerance on the task regulation: a task is considered as regulated when e i (t) A task sequence is characterized by a set of timeconstraints binding the schedules of two tasks e i and e j . They can be defined as follow2 : e i must begin or end once e j has begun, has ended or has been regulated.

≤ ǫ i . The regulation time t R i is defined by e i (t R i ) = ǫ i . j begins once i has begun t I i ≤ t I j j begins once i is realized t R i ≤ t I j j begins once i has ended t F i ≤ t I j j ends once i is realized t R i ≤ t F j j ends once i has ended t F i ≤ t F j
We use the graphical representation given by Fig. 1 and the following notation to describe the sets of pairs of tasks e i and e j that undergo these dependencies (e i is the direct predecessor of e j ) :

S I,I = {(e i , e j ) | t I i ≤ t I j } (5a) S R,I = {(e i , e j ) | t R i ≤ t I j } (5b) S F,I = {(e i , e j ) | t F i ≤ t I j } (5c) S R,F = {(e i , e j ) | t R i ≤ t F j } (5d) S F,F = {(e i , e j ) | t F i ≤ t F j } (5e)
For example, the robot has first to grasp an object and maintain the force closure on it (e A ) before moving it (e B ). The task (e B ) can only start once the task (e A ) has been realized, and must end before the task (e A ).

III. CONTINUOUS OPTIMIZATION OF SEQUENCE OF

TASKS

Given a set of hypothesis described using [START_REF] Ghallab | Automated Planning: Theory and Practice[END_REF], we now propose a generic solution to automatically compute an optimal set of task-behavior parameters and their sequencing plan to be executed by the SoT.

A. General problem formulation

An optimization problem is composed of a criterion to minimize, and of a set of equality and inequality constraints that must be satisfied. Our chosen criterion is to minimize the regulation duration of the mission. The variables of our problem are for each task: (i) the time of its entry, (ii) the time of its removal (from the SoT), and (iii) the gains (λ I , λ F , β) which describe the task execution behavior.

The general optimization problem is written as follows:

min x t End (6a) subject to q = SoT x (q, t) (6b) seq(q) < 0 (6c) φ(q) < 0 (6d) ∀i, t F i ≤ t End (6e) The vector x = [t I 1 , t F 1 , λ I 1 , λ F 1 , β 1 , . . . , t I n , t F n , λ I n , λ F n , β n , t
End ] gathers the optimization variables of each task and t End , the duration of the mission. seq(q) and φ(q) are respectively the sequencing and the robotic constraints.

The optimization criterion t End is computed indirectly. An equivalent explicit definition could be given by t End = max i (t F i ). However this constraint is not smooth. Giving only (6b), the problem is smooth and properly defined: at the optimal solution, t End will be equal to the maximum termination time of all tasks' t F i . Vector q is in fact a vector of functions of time, hence constraints φ(q) are semi-infinite, i.e. taking place for all the values of the continuous variable t ∈ [t 0 , t End ].

It can be shown that ( 6) defines a continuous optimization problem. However, it cannot be solved directly because of the semi-infinite nature of the constraints. Therefore we expanded the semi-infinite constraint into a discreet form.

B. Constraints

Parameter x must satisfy both the sequencing and the robotic time-constraints enumerated hereafter:

1) Tasks constraints, noted seq(q): gather the task sequence conditions of (5) and the following constraints: For each task i:

Time coherence 0 ≤ t I i < t F i ≤ t End (7a) Termination condition s * i -s i (t F i ) < ǫ i (7b) Gain consistency λ I i ≤ λ F i (7c)
The constraints (5a), (5c), (5e), (7a) and (7c) are linear. On the contrary, the constraint (7b) is impossible to compute directly using x, and is determined from a simulation of the execution. Care has to be taken while resolving the condition described by (5b) and (5d). Indeed, discretizing t R to the closest simulation step will produce discontinuities which may disturb the optimization process. A rather fastidious solution to this continuity problem would be to determine this point by interpolation. Another solution is to reformulate them by evaluating the regulation of the task i instead. The constraint (5b) and (5d) becomes respectively:

∀(i, j) ∈ S R,I , s * i -s i (t I j ) ≤ ǫ i (8a) ∀(i, j) ∈ S R,F , s * i -s i (t F j ) ≤ ǫ i (8b)
2) Robot constraints : φ(q): Those constraints are mainly due to hardware intrinsic limitations of the robot: Joint limits q min ≤ q ≤ q max (9a) Velocity limits qmin ≤ q ≤ qmax (9b)

Collision avoidance 0 ≤ d ij (9c)
q min , q max , qmin , qmax are respectively the lower and upper joint limits and the lower and upper velocity limits. d ij is the distance between objects i and j. Object designate those found in the mission's environments and each link of the robotic system. Hence, both collision with the environment and self-collision of the robot have to be evaluated.

All of those constraints are semi-infinite: the following section presents how they have been tackled.

C. Technical aspects of the optimization resolution 1) Semi-infinite constraints:

In a first approach, we tried to discretize the semi-infinite constraints on the basis of the simulation steps grid. However, since the number of the grid sample points changes in function of t End , the number of constraints is variable. Subsequently a classical optimization solver can not handle them.

Let c be the evaluation value of a given constraint: (∀t ∈ [t 0 , t End ], c(t) < 0). We considered associating only one value to the constraint, c V , that is computed as follows: If the constraint is always satisfied, then c V is the higher value of c(t). Otherwise, it is the sum of all the violations found at each time step. Considering that the time step can change (e.g. when adding an interpolation point), we choose to weight the added value by the time step δt.

2) Constraint by task: Each task appears only once in the sequence, but a same action can be associated to many tasks. Associating the constraints φ(q) to the whole simulation can thus raise an issue: a violated constraint can not be linked to the responsible task. In order to compensate this problem, we consider n T additional sets of constraint φ(q), noted φ i (q), i ∈ [1 . . . n T ], (with n T the number of tasks in the sequence). Each set φ i (q) is computed only when the task i is in the SoT.

3) Scaling: Since the constraints are not homogeneous (times, angles, velocities, distances), they have to be normalized based on the constraint values obtained while executing the sequence corresponding to the initial set of parameters x 0 . This simple scaling improves significantly the convergence of the optimization.

D. Absolute versus relative timing

In this parameterization, the tasks are described with an absolute time. As it is, decreasing t I i for a task i will not have any direct effect on t F i : we have also to decrease t F i then decrease t End : it is thus necessary to propagate the reduction for all the following tasks. To avoid this, another parameterization consists in describing the SoT entry time of a given task with respect (i.e. relatively) to the previous one.

We introduce a relative timing: each task is now described by two delays (instead of the absolute times t I and t F ), namely: 1) dt I : is the delay which occurs between (i) the maximum time of entry or of end of the preceding tasks, and (ii) the SoT entry time of the task in question.

t I i = max max (j,i)∈SI,I {t I j }, max (j,i)∈SF,I {t F j } + dt I (10)
2) dt F : is the delay between the SoT entry and the removal times of the task in question.

t F i = t I i + dt F i (11)
Subsequently, the new parameter vector is noted :

x ′ = [dt I 1 , dt F 1 , λ I 1 , λ F 1 , β 1 , . . . , dt I n , dt F n , λ I n , λ F n , β n , t End ].
If the task sequence is only a chain of tasks realized one after the other, we directly have x ′ = f (x), with f a linear function, and

t End = i (dt I i + dt F i )
Considering this new set of parameters, the formulation of the optimization problem changes: some tasks constraints of seq(q) are modified. The previous constraints (5a), (5c) and (7a) are replaced by these constraints on the delay:

∀i,0 ≤ dt I i (12a) ∀i,0 < dt F i (12b) IV. IMPLEMENTATION

A. Optimization

At each optimization step, the solver chooses a new set of parameters x. It then computes the constraints. Constraints (5e) and (7c), (12a) and (12b), can be evaluated directly. As stated previously, the other constraints can not be directly computed (since they do not write in an analytical formulation). They are thus evaluated using a complete simulation of their execution. The chosen value of the current optimization variable vector x is transmitted by the optimization solver to the simulation engine. The simulation returns the evaluation of the constraints and the optimization solver computes a new step vector x, until convergence.

Our optimization problem is a non-linear contrained parametric problem. We chose the SQP algorithm from the MATLAB optimization toolbox, which is suitable to this kind of problem.

B. Simulation

In section II, we presented the computation of desired joint velocities for a hierarchy of tasks, as [START_REF] Hanafusa | Analysis and control of articulated robot with redundancy[END_REF]. The simulation is basically a numerical integration of this equation (we used an explicit Euler integration method with a fixed step ∆t = 0.005sec). The entry and exit times t I i and t F i are continuous variables that are not aligned with the grid. Those instants are important since they correspond to a change in the SoT and thus a change in the control. If postponing the change of control to the next time step (like on a real system) we will not have a continuous problem (hence potentially raising the same problem described in section III-B). The algorithm 1 describes the simulation. The function computeTimes computes the absolute times using the relative times. The function findTimeStep computes the required time step for the Euler integration: the initial ∆t, or a smaller one if needed, due to the need of splitting this interval in two. The function handleStackOfTasks computes the velocity of the robot induced by the tasks execution and integrates it, altogether with any other simulated objects or processes, to obtain the new positions.

The simulation engine runs under the AMELIF framework [START_REF] Evrard | Framework for haptic interaction with virtual avatars[END_REF], an interactive dynamic simulator for virtual avatars which includes collision detection and task handling according to the SoT formalism. The execution for both simulation and real-robot control is performed by a generic control framework based on [START_REF] Mansard | Task sequencing for sensor-based control[END_REF].

V. EXPERIMENT

A. Temporal network

The sequence of tasks (Fig. 2) describes a robot taking out a can from the fridge. The corresponding tasks are: The task e 6 is an intermediary task introduced as a way point: its tolerance on task regulation ǫ 6 is large so that the arm does not have to stop. This is part of the optimization decision, in order to reduce the execution time. This is a complex mission that can not be split into smaller sequences. Indeed, the sequence is centered on the fridge: the grasping part does not make sense if the fridge is closed. Instead of adding explicit timing conditions between the tasks to ensure that this will never occur, we choose to consider as constraint the collision between the left arm and the door, in order to allow task overlapping.

The constraints considered for this problem are thus sequencing and robotic constraints (joint position and velocity limits), and collision avoidance with the fridge. Fig. 2. Sequence describing the HRP-2 taking the can in the fridge Fig. 3. Results of the optimization of the sequence of task: when the task is added in the SoT, its error is first regulated (this is the dark part (red or dark blue) of the block). From t R i , the error is nearly null and the task is kept in the SoT (light part (yellow or cyan) of the block) until t F i .

B. Results of the optimization

We ran the optimization on a 3GHz desktop PC running under Windows OS. The sequence found is described on Fig. 3. Each task is described by two periods: the dark one is the achievement period [t I i , t R i ], the bright one is the SoT presence period [t I i , t F i ]. The overlaps between the tasks of the left and the right arm appear clearly: the left arm starts to move before the fridge is open. It then starts to move toward the can pose even if the fridge is not completely open. And finally, the right arm starts to close the fridge before the left arm has completely left the fridge area. The whole task sequence lasts 47sec. Without these two overlaps, the robot will move to and grasp the can (e 7 ) only after the fridge is fully opened (e 3 ) and it will close the fridge (e 4 ) only after the can is completely taken out (e 10 ); consequently the total mission would have taken at least 71sec.

C. Experiment on the real robot

The task sequence is experimented on the upper body of the HRP-2 humanoid robot. Only the described tasks are used to compute the control law (which means that no additional care is taken for ensuring the constraints). For the tasks requiring a haptic interaction (i.e. opening and closing the fridge) the force sensor of the robot is used to close the loop and compensate for position uncertainties.

The robot manages to grasp the can without colliding any obstacle or joint limits, and respecting the given velocity limits. The obtained execution is plotted on Fig. 4. Thanks to the optimized gain, the convergence of the error of the tasks that require a good precision (grasping the fridge handle and the can) is achieved as quickly as allowed by joint velocity limits. Snapshots of the execution are given in Fig. 5

VI. CONCLUSION

We devise a method which allows to optimize both the behavior and the overlapping scheduling of a sequence of tasks composing a robotic mission. The solution derives from an optimization formulation of the tasks scheduling keeping the formalism built on the top of a task-function based control. This allows to include the robot limitations as well as collision avoidance as constraints. Our method is exemplified through a complete simulation of a complex mission, where we demonstrated an improvement in the smoothness of the generated motion. For the time being, our method still needs a predefined ordered sequence. As a future work we will increase the autonomy by determining automatically the

Fig. 1 .

 1 Fig. 1. Five time-dependency relations are considered.

•

  Tasks of the right arm: (e 0 ): open the gripper, (e 1 ): move it to the fridge handle, (e 2 ): close it, (e 3 ): open the fridge, (e 4 ): close the fridge • Tasks of the left arm: (e 5 ): open the gripper, (e 6 ): move it in the fridge area, (e 7 ): move it to the can, (e 8 ): close it, (e 9 ): lift the can, (e 10 ): remove the can out of the fridge,

Fig. 4 .

 4 Fig. 4.Experiment on HRP-2: errors diminish when optimized task scheduling is applied: (top) right arm tasks (middle) left arm tasks (bottom) gripper tasks. The concurrency between the tasks is clearly visible.

  3 .

Fig. 5 .

 5 Fig. 5. HRP-2 grasping a can in the fridge.

  To solve this problem, the entry time t a is added as an integration point during the time step [t, t + ∆t], splitting it into the two smaller ones [t, t a ] and [t a , t + ∆t].

	Initialization
	[t I 1 , t F 1 , . . . , t I n , t F n ] = computeTimes (x) t End Sim = max i t F i , t = 0
	while (t < max (t End , t End Sim )) do ∆t ′ = findTimeStep(t)
	handleStackOfTasks (t)
	updateConstraints()
	t = t + ∆t ′ end
	Algorithm 1: Tasks sequencing simulation

contrary to Allen Logic, that only considers the start and end points of the time interval, here is also considered the regulation time t R

www.laas.fr/ ˜fkeith/iros09.avi, the video is also attached to the paper ordered sequence and compute all the necessary subtasks from definitions of actions/objects associations. We will also focus on more complex scenario using in particular perception tasks such as visual interaction.

ACKNOWLEDGMENT

This work is partly supported by the European project FP6 ROBOT@CWE www.robot-at-cwe.eu.