A P Rinaldi
email: aprinaldi@lbl.gov

ISTerre J Vandemeulebrouck
email: jean.vandemeulebrouck@univ-savoie.fr

M Todesco
email: micol.todesco@bo.ingv.it

F Viveiros

Effects of atmospheric conditions on surface diffuse 1 degassing 2

Diffuse degassing and fumarolic activity are secondary manifestations of active volcanoes. Monitoring of these manifestations plays a key role in understanding the activity status of volcanoes, and the installation of permanent stations to measure the CO 2 flux at the surface are becoming more common. The measured gas fluxes often present a correlation with environmental observables, such as air temperature or barometric pressure, especially when gas emission is low. An example is the Furnas volcano (Azores), which features low temperature fumaroles, CO 2 -rich springs, and several diffuse degassing areas. The hydrothermal CO 2 flux through the soil is continuously recorded by permanent gas stations coupled with several meteorological sensors. Daily and seasonal cycles have been observed in the time series of CO 2 flux. Air temperature and barometric pressure are the variables Vulcanologia, Sezione di Bologna, Italy.

that best correlate with the CO 2 flux cycles. To understand the influence of atmospheric and soil conditions on the gas emission, several simulations with TOUGH2 simulator were performed. A simple 1D model was developed to understand the physical mechanisms producing the observed variations. Then, a dual parameters study focused on how rock permeability and gas source properties affect the resulting fluxes. Numerical results are in agreement with the observed data and show that changes in the barometric pressure and air temperature may cause changes in the observed CO 2 fluxes, which depend on rock permeability and on pressure driving the upward flux.

Introduction

Continuous measurements of carbon dioxide flux are now commonly used to monitor the degassing in volcanic environments [START_REF] Chiodini | Soil CO 2 flux measurements in volcanic and geothermal areas[END_REF][START_REF] Chiodini | Long term variations of the Campi Flegrei (Italy) volcanic system as revealed by the monitoring of hydrothermal activity[END_REF][START_REF] Hernández | Carbon Dioxide Degassing by Advective Flow from Usu Volcano, Japan[END_REF].

The transfer of the CO 2 to the surface results from two different physical processes, namely advection and diffusion [e.g. [START_REF] Werner | CO 2 emission related to the Yellowstone volcanic system 2. Statistical sampling, total degassing, and transport mechanisms[END_REF]. Advective transport is mainly induced by pressure gradient, and by temperature gradients that create density differences. Diffusion, on the other hand, is controlled by concentration gradients. For both advection and diffusion processes, the flux depends on the properties of the gas source, on the petrophysical properties of the soil, e.g. permeability, and on the conditions at the surface. In this paper, we aim at modeling how the periodic variations of air temperature and pressure at the soil surface can affect the CO 2 flux signals by producing harmonic variations at the same periods. A second goal is to understand if the monitoring of the periodic flux variations could be used to detect source or medium properties changes. Our model will be tested on CO 2 flux signal recorded at Furnas volcano (Azores, Portugal).

The impact of environmental variables on soil gas emissions has been studied in several degassing areas in the last decades [e.g. [START_REF] Clements | Atmospheric pressure effects on 222 Rn transport across the Earth-Air interface[END_REF][START_REF] Reimer | Use of Soil-Gas Helium Concentrations for Earthquake prediction: Limitations Imposed by Diurnal Variation[END_REF][START_REF] D R A F T Klusman | Preliminary analysis of meteorological and seasonal influences on crustal gas emission relevant to earth prediction[END_REF][START_REF] Hinkle | Factors affecting concentrations of helium and carbon dioxide in soil gases[END_REF][START_REF] Hinkle | Environmental conditions affecting He, CO 2 , O 2 , and N 2 in soil gases[END_REF][START_REF] Chiodini | Soil CO 2 flux measurements in volcanic and geothermal areas[END_REF][START_REF] Mcgee | Annual cycle of magmatic CO 2 in a tree-kill soil at Mammoth Mountain, California: Implications for soil acidification[END_REF][START_REF] Mcgee | Geochemical evidence for a magmatic CO 2 degassing event at Mommoth Mountain, California[END_REF][START_REF] Rogie | Dynamics of carbon dioxide emission at Mammoth Mountain, California[END_REF][START_REF] D R A F T Granieri | Continuous monitoring of CO 2 soil diffuse degassing at Phlegrean Fields (Italy): influence of environmental and volcanic parameters[END_REF][START_REF] Granieri | Carbon dioxide diffuse emission from the soil: ten years of observations at Vesuvio and Campi Flegrei (Pozzuoli), and linkages with volcanic activity[END_REF][START_REF] Hernández | Radon and helium in soil gases at Canãdas Caldera, Tenerife, Canary Islands, Spain, J. Volcanol[END_REF][START_REF] Lewicki | Dynamic coupling of volcanic CO 2 flow and wind at the Horseshoe Lake tree kill, Mammoth Mountain, California[END_REF][START_REF] Viveiros | Environmental influences on soil CO 2 degassing at Furnas and Fogo volcanoes (São Miguel Island, Azores archipelago)[END_REF][START_REF] Viveiros | Meteorological factors controlling soil gases and indoor CO 2 concentration: a permanent risk in degassing areas[END_REF][START_REF] Cigolini | Radon surveys and real-time monitoring at Stromboli volcano: Influence of soil temperature, atmospheric pressure and tidal forces on 222 Rn degassing[END_REF]. However very few papers were found in the literature that analyzed diurnal variations for the soil CO 2 flux in hydrothermal environments. [START_REF] Chiodini | Soil CO 2 flux measurements in volcanic and geothermal areas[END_REF] noticed that CO 2 flux variations at Vulcano (Aeolian Islands, Italy) were affected by barometric pressure changes D R A F T August 30, 2012, 9:29am D R A F T rather than by absolute atmospheric pressure value. [START_REF] D R A F T Granieri | Continuous monitoring of CO 2 soil diffuse degassing at Phlegrean Fields (Italy): influence of environmental and volcanic parameters[END_REF] identified a 24 h cycle in CO 2 fluxes from Solfatara volcano between 1998[START_REF] Aumento | Radon tides on an active volcanic island: Terceira, Azores[END_REF][START_REF] Padrón | Diffuse CO 2 emission rate from Pululahua and the lake-filled Cuicocha calderas, Ecuador[END_REF] recognized diurnal and semidiurnal fluctuations in the gas flux acquired during 2004 in a station located at El Hierro (Canary Islands). On the other side, diurnal and seasonal variations of CO 2 flux derived from soil respiration (biogenic origin) have been intensely studied as CO 2 production is dependent on the temperature and the CO 2 oscillations are positively correlated with the air/soil temperature changes [START_REF] Witkamp | Environmental effects on microbial turnover of some mineral elements: Part II-Biotic factors[END_REF][START_REF] Bajracharya | Diurnal and seasonal CO 2 -C flux from soil as related to erosion phases in central Ohio[END_REF][START_REF] D R A F T Nakadai | Diurnal changes of carbon dioxide flux from bare soil in agricultural field in Japan[END_REF], as well as with the wind speed [START_REF] Takle | Influence of high-frequency ambient pressure pumping on carbon dioxide efflux from soil[END_REF][START_REF] Reicosky | Tillage and wind effects on soil CO 2 concentrations in muck soils[END_REF][START_REF] Bowling | Persistent wind-induced enhancement of diffuse CO 2 transport in a mountain forest snowpeak[END_REF]. Biogenic fluxes are also influenced by atmospheric pressure [START_REF] Massman | Advective transport of CO 2 in permeable media induced by atmospheric pressure fluctuations: 1. An analytical model[END_REF]. These papers recorded the correlation but did not investigate the processes that cause and explain the observed trends.

Some studies have also recognized periodicities on 222 Rn time series in the air [e.g. [START_REF] Rigby | Wind and barometric pressure effects on Radon in two mitigated houses[END_REF][START_REF] Pinault | Signal processing of diurnal and semidiurnal variations in radon and atmospheric pressure: A new tool for accurate in situ measurement of soil gas velocity, pressure gradient, and tourtuosity[END_REF][START_REF] Robinson | Soil-gas entry into houses driven by atmospheric pressure fluctuations -The influence of soil properties[END_REF][START_REF] Groves-Kirkby | Identification of tidal and climatic influences within domestic radon time-series from Northamptonshire, UK[END_REF][START_REF] Steinitz | Annual cycle of magmatic CO 2 in a tree-kill soil at Mammoth Mountain, California: Implications for soil acidification[END_REF][START_REF] Richon | Detectability and significance of 12 hr barometric tide in radon-222 signal, dripwater flow rate, air temperature and carbon dioxide concentration in an underground tunnel[END_REF] and in the soil [START_REF] Aumento | Radon tides on an active volcanic island: Terceira, Azores[END_REF][START_REF] Richon | Radon anomaly in the soil of Taal volcano, the Philippines: A likely precursor of the M7.1 Mindoro earthquake[END_REF][START_REF] Richon | Temporal signatures of advective versus diffusive radon transport at a geothermal zone in Central Nepal[END_REF][START_REF] Cigolini | Radon surveys and real-time monitoring at Stromboli volcano: Influence of soil temperature, atmospheric pressure and tidal forces on 222 Rn degassing[END_REF]. For example, [START_REF] Clements | Atmospheric pressure effects on 222 Rn transport across the Earth-Air interface[END_REF] observed that pressure variations of 1% -2% changes the radon flux from 20% to 60%.

This work intends to understand the influence of atmospheric and soil conditions on the gas release, by performing several simulations with TOUGH2 geothermal simulator [START_REF] Pruess | TOUGH2 user's guide, version 2.0[END_REF]. Although such a simulator is largely used in non-volcanological contests [e.g. [START_REF] Oldenburg | Buoyancy Effects on Upward Brine Displacement Caused by CO 2 Injection[END_REF][START_REF] Borgia | Numerical simulation of salt precipitation in the fractures of a CO 2 -enhanced geothermal system[END_REF][START_REF] Mazzoldi | Induced seismicity within GCS projects: maximum earthquake magnitude and leakage potential from undetected faults[END_REF],

applications in volcanology are becoming more common for the study of diffuse degassing [e.g. [START_REF] Todesco | Modeling of unrest signals in heterogeneous hydrothermal systems[END_REF][START_REF] Chiodini | Long term variations of the Campi Flegrei (Italy) volcanic system as revealed by the monitoring of hydrothermal activity[END_REF] and of the geophysical and geochemi-

D R A F T August 30, 2012, 9:29am D R A F T
cal signals related to hydrothermal fluids circulation [e.g. [START_REF] Hutnak | Numerical models of caldera deformation: Effects of multiphase and multicomponent hydrothermal fluid flow[END_REF][START_REF] Rinaldi | Hydrothermal instability and ground displacement at the Campi Flegrei caldera[END_REF][START_REF] Rinaldi | Electrical conductivity, ground displacement, gravity changes, and gas flow at Solfatara crater (Campi Flegrei caldera, Italy): results from numerical modeling[END_REF]. We used the TOUGH2/EOS2 module to describe CO 2 in gas phase fluid. Using a 1D model, a parametric study was performed to understand the physical mechanisms producing the observed variations. Numerical results, in agreement with the observed data, show that the CO 2 fluxes are strongly dependent on reservoir pressure, temperature and pressure changes applied at the surface, and on domain permeability. with the periodical behavior of some meteorological parameters. Statistical analysis applied to the recorded time series shows that air temperature, barometric pressure, rainfall, and wind speed are the observables that better correlate with the soil CO 2 flux variations [START_REF] Viveiros | Environmental influences on soil CO 2 degassing at Furnas and Fogo volcanoes (São Miguel Island, Azores archipelago)[END_REF][START_REF] Viveiros | Meteorological factors controlling soil gases and indoor CO 2 concentration: a permanent risk in degassing areas[END_REF]. These external parameters may also have a different effect on the gas flux depending on the characteristics of the monitoring site (e.g. physical properties of the soil, topographic effects, drainage area). [2010], the stations are all installed in diffuse degassing structures (DDS) that are fed by the hydrothermal sources.

A case

The permanent flux stations perform measurements based on the accumulation chamber method [START_REF] Chiodini | Soil CO 2 flux measurements in volcanic and geothermal areas[END_REF]. Every hour, a chamber is lowered on the ground and the gas released at the ground surface is pumped into an infrared gas analyzer (IRGA).

The soil CO 2 flux value is computed as the linear best fit of the flux curve during a predefined period of time. Measurements have a reproducibility within 10% for the CO 2 range between 10 and 20000 g m -2 d -1 [START_REF] Chiodini | Soil CO 2 flux measurements in volcanic and geothermal areas[END_REF]]. The gas flux automatic stations have also coupled meteorological sensors and acquire simultaneously the barometric pressure, air temperature, air relative humidity, wind speed and direction, rainfall, soil water content and soil temperature.

2.1. Spectral analysis applied to CO 2 flux data at GFUR2

Seasonal periodicities and diurnal variations are observed on the CO 2 time series concomitant with the periodicity observed in some meteorological variables (e.g. air temperature, air relative humidity, and barometric pressure). We choose to analyze data

D R A F T August 30, 2012, 9:29am D R A F T
registered at station GFUR2 during the summers of the whole period (2005 -2009, from May to September), since, in that period of the year, average monthly rainfall is low, and then rainfall effects on degassing can be considered as negligible. In addition, GFUR2 is the monitoring site with longer register at Furnas volcano.

The magnitude of the diurnal variations on the soil CO 2 flux (usually from 50 to 100 g m -2 d -1 ) excludes the possibility of a biogenic origin of the CO 2 variations. Works published in biogenic environments [e.g. [START_REF] D R A F T Nakadai | Diurnal changes of carbon dioxide flux from bare soil in agricultural field in Japan[END_REF] refer significantly lower daily amplitudes (about 5 -6 g m -2 d -1 ), indicating that the site of the present study at Furnas are clearly fed by volcanic-hydrothermal sources.

Diurnal and semi-diurnal components are observed in the soil CO 2 flux time series (Fig.

2

), with 24h period peak significantly stronger than 12h period peak. Air temperature, air relative humidity, barometric pressure, and wind speed recorded at GFUR2 site show as well daily fluctuations. Figure 3 shows an example of a summer week variation of CO 2 flux, air temperature, and barometric pressure, observed at GFUR2. The diffuse emission of carbon dioxide reaches maximum values early in the morning, when atmospheric pressure and temperature are lower. The variations of the two atmospheric parameters will be used as an evolving, top boundary condition in the numerical simulation (see next section).

The correlation coefficient between the monitored variables and the time delays between them were calculated, for the 24h and 12h peaks. Only correlation coefficients higher than 0.5 were included in Table 1. For the 24 h component, the soil CO 2 presents the highest correlation (r > 0.8) with wind speed, air temperature, and air relative humidity. With respect to the 12h component, air temperature and wind speed are the meteorological variables that better correlate with the CO 2 flux (r ∼ 0.7). We considered the wind as a D R A F T August 30, 2012, 9:29am D R A F T perturbation of lower period (less than 12 h), which means that its effect is very shallow, with a penetration up to some decimeters but with a large impact. [START_REF] Viveiros | Environmental influences on soil CO 2 degassing at Furnas and Fogo volcanoes (São Miguel Island, Azores archipelago)[END_REF] recognized the inverse correlation between CO 2 flux and wind speed and suggested that during high wind speed events the soil gas is diluted with atmospheric air pushed into the upper parts of the soil. This correlation is potentially favored by the soil porosity and permeability. Barometric pressure shows in general a worse correlation with diffuse degassing, with the lowest coefficient found for the 24h component, but it better correlates with the soil degassing for the 12h component (r ∼ 0.7). In our case, the CO 2 flux released by an hydrothermal area has an opposite behavior to that observed where CO 2 is entirely biogenic, and where the CO 2 oscillations are positively correlated with temperature changes [START_REF] Witkamp | Environmental effects on microbial turnover of some mineral elements: Part II-Biotic factors[END_REF][START_REF] Bajracharya | Diurnal and seasonal CO 2 -C flux from soil as related to erosion phases in central Ohio[END_REF][START_REF] D R A F T Nakadai | Diurnal changes of carbon dioxide flux from bare soil in agricultural field in Japan[END_REF].

Numerical simulation

In order to compute the CO 2 flux changes that arise from variations in atmospheric pressure and temperature conditions and to understand the role played by rock permeability and reservoir overpressure on observed delay, parameter studies of fluid circulation were carried out with the TOUGH2/EOS2 simulator [START_REF] Pruess | TOUGH2 user's guide, version 2.0[END_REF]. Briefly, TOUGH2 is

D R A F T August 30, 2012, 9:29am D R A F T
a multi-purpose simulator for fluids in porous medium. Using an integral finite difference method for mass and energy balance, and a first order finite difference for time discretization, fluids can be simulated as multi-phase and multi-component, taking into account both the effects of relative permeability of each phase and capillarity pressure. The module EOS2, here used to simulate the CO 2 flux only, may account for the presence of water and carbon dioxide. Heat transfer in TOUGH2 occurs both by conduction and convection, and it accounts for latent heat effects. A thermodynamic equilibrium is assumed to be present between fluids and porous matrix. We do not account for chemical reaction nor deformation of the rock matrix. At this time, we focus on the unsaturated, upper portion of the soil that we consider fully saturated with carbon dioxide. We are aware that the presence of liquid water may significantly alter the simulated CO 2 flux, but this effect is beyond the limits of the present work. The equations solved in TOUGH2/EOS2

for single-phase conditions and CO 2 as the only fluid component are shown in Table 2.

We aim to analyze the effects of the atmospheric conditions on diffuse degassing looking at diurnal and semi-diurnal variation of CO 2 flux, and for this reason we simulated a weeklong variation of the air temperature and barometric pressure, modifying the boundary conditions at the top of the domain. The variations of the physical quantities correspond to data collected at the GFUR2 station (Fig. 3). Initial conditions correspond to a steady state that was reached with a long simulation (hundreds of year), which produced a linear pressure and temperature gradients in the medium.

Furnas case

We first present a detailed analysis of a base case, in order to describe which mechanisms drive and affect the degassing. Rock properties and bottom boundaries conditions for this base case were chosen to match the CO 2 flow variation observed at GFUR2 during the considered week (Fig. 3). We chose a domain permeability κ = 2•10 -14 m 2 and an overpressure in the reservoir ∆p = 0.05 MPa, with respect to the initial atmospheric pressure. At steady state conditions, the system is characterized by temperature and pore pressure linear gradients, and by a stationary upward CO 2 flow. Starting from these conditions, we run a week-long simulation, imposing the observed variations of atmospheric pressure and temperature along the top boundary. Results of modeled CO 2 flux at the surface obtained with the chosen permeability and gas reservoir overpressure are in good agreement with the changes observed at GFUR2 (Fig. 5b). Some discrepancy may result from the fact that we do not take into account effects of wind, air humidity, and soil water content.

The effects of the evolving boundary conditions are not confined to the ground surface, but propagate downward, causing changes of the CO 2 flux throughout the domain (Fig. 5c). The CO 2 flux observed at the base of the domain presents almost the same variation of the degassing simulated at ground surface.

The applied pressure and temperature perturbations propagate over the 1 m-high domain (Fig. 6a,b,c,d). Applied atmospheric pressure changes are of the order of few hundreds Pa, and this pressure perturbation often reaches most of the computational domain (Fig. 6b). However, these changes are small compared to the average pressure within the system (order of 10 5 Pa). Although the initial pressure distribution is mostly unaffected (Fig. 6a), the changes in the atmospheric pressure clearly influence the simulated CO 2 flux. On the other hand, the changes of air temperature (Fig. 6d) are larger compared

D R A F T August 30, 2012, 9:29am D R A F T
to the average temperature of the whole system, and clearly alter the initial temperature distribution (Fig. 6c). Both pressure and temperature perturbations cause changes of gas density (Fig. 6e,f).

Comparison between harmonic flux changes in model and observations

Power spectra

In order to understand which atmospheric perturbation has the largest effect on the resulting CO 2 degassing, a spectral analysis has been performed. When a perturbation of a certain period generates large variations in the resulting time series, then the data will show a large peak in the spectrum at that particular period.

Using as example the Furnas case analyzed in the previous section, we now analyze the power spectrum of the gas flux time series to find out its characteristic frequencies.

In order to focus on the diurnal and semi-diurnal components, we low-pass filtered the data of the week-long record, and keep only periods lower than 28 h (Fig. 7). Without longer periods, the pressure has a main 12h period, with a very small component at 24h.

The temperature is only slightly affected by the filter, since it is mainly characterized by periods shorter than 24h. Temperature has a main peak at 24h, and a twenty-times smaller peak at 12h. The gas flux in this base case has a main 24h period, but the 12h component is strong enough to affect the temporal evolution. The ratio 12h peak/24h peak for the gas flux is about 7%, larger than the 5% we found for the temperature. This suggests that in this base case, the temperature variation mostly affect the evolution of gas flux (for the 24h period). A small variation in the 12h period is noticeable, and it is due to the application of the pressure perturbation, otherwise the ratio 12h peak/24h peak would have been the same for gas flux and temperature. The influence of pressure

D R A F T August 30, 2012, 9:29am D R A F T
can also be seen in the effect of long-term (several days) pressure changes, as during the studied week where pressure and CO 2 long-term changes are well correlated.

Time delay

The time delay between the applied perturbations (atmospheric pressure or temperature) and the resulting CO 2 flux time series may provide interesting information on the properties of the soil and on the conditions of the hydrothermal reservoir. We apply a band pass filter centered on the selected period (24h or 12h) with a 4 h bandwidth, in order to compare only the 24h component of the gas flux to the air temperature, and only the 12h component of the gas flux to the atmospheric pressure. It is of note that the 12h peak of the power spectrum is one order of magnitude below the 24h.

Results of the band-pass filtering of the simulated gas flux are shown in Figure 8 and compared with the filtered pressure (Fig. 8b) and temperature (Fig. 8c), respectively.

The time series are now very smooth and comparable, and it is quite easy to calculate a time delay as a time shift between the two series by a correlation method. For the base case simulation (normalized resulting gas flux and perturbations in Fig. 8a), results

show a time delay of 6 hours when compared with the atmospheric pressure time series (Fig. 8b), and a time delay of 15.25 hours when compared with the temperature time series (Fig. 8c). Time delay may change when domain properties are changed (see next section), however, these results are in good agreement with the delays observed at Furnas The first set of simulations was run changing the rock permeability and keeping a constant gas reservoir overpressure as in the base case (0.05 MPa, a value that is consistent with observed gas flux). The resulting changes in peak amplitude and time delay are shown in Figure 9. The peak amplitude of the simulated CO 2 flux always presents both the 12h and 24h component (Fig. 9a, green and blue lines, respectively). However, the two analyzed components are not always comparable. The 24h component is dominant for most of the considered permeabilities. More than one order of magnitude differences between 24h and 12h components appear for low permeabilities. This difference decreases at higher permeability, and the two considered components have the same order of magnitude for the highest permeability considered (10 -11 m 2 ). Results also show changes in time delay (Fig. 9b). The 24h component of the CO 2 flux has a delay with respect the temperature perturbation ranging between 18 hours and 14 hours (Fig. 9b, blue line).

The higher the permeability, the lower the delay observed. The delay of the 12 h component of the CO 2 flux with respect to the pressure perturbation is rather constant (∼ 6 hours) for most of the simulated permeabilities, and increases up to 8 hours only for permeability lower than 10 -14 m 2 (Fig. 9b, green line).

D R A F T August 30, 2012, 9:29am D R A F T
The second set of simulations was run keeping a constant rock permeability, set as in the base case (2•10 -14 m 2 ), and changing the overpressure at the base of the domain. Figure 10 shows the changes in peak amplitude and time delay as a function of the gas reservoir overpressure. Again, all the simulated cases present both the 12h and 24h components.

The peak amplitude for the 12 h component of the CO 2 flux has a constant value for low overpressure (up to 0.01 MPa), then it exponentially decreases for higher overpressure (Fig. 10a, green line). The 24h component features an opposite behavior: it starts to increase exponentially for overpressure greater than few thousands Pa (Fig. 10a, blue line). The ratio between the peak amplitudes of the two considered components reaches a maximum of three orders of magnitude for the highest simulated overpressure. This means that for low overpressure the system seems to have only a 12h component, which indicates that it is only affected by pressure perturbations. At higher overpressures, the CO 2 flux appears to have only a 24h component, which means that only the temperature perturbation affects the resulting degassing. The time delay only slightly changes as a function of the gas reservoir overpressure. The 12h component has a constant delay of ∼6 hours with respect to the pressure for all the considered overpressure (Fig. 10b, green line). The time delay of the 24h component with respect to the temperature is larger than 18 hours for low overpressure, and starts to decrease when the overpressure is greater than one thousand Pa, up to a minimum, constant value of ∼15 hours for the highest overpressures simulated (0.05 -0.1 MPa) (Fig. 10b, blue line).

Figure 11 shows a dual parameters analysis of the time delay and peak amplitude in the spectrum. Throughout the studied parameters domain, the 12h component of the CO 2 flux presents a stable delay around 6 hours with respect the barometric pressure (Fig.

D R A F T August 30, 2012, 9:29am D R A F T X -17
11a), and reaches a maximum of 8 hours only for very low permeability (10 -15 m 2 ). The delay of the 24h component with respect the atmospheric temperature is shown in Figure 11b. This delay remarkably changes in the parameters space: a high delay (∼ 21 hours) is observed for low permeability and low gas reservoir overpressure. Then, this time delay lowers when both permeability and gas reservoir overpressure increase, the minimum value (∼ 13 hours) being observed for the highest gas reservoir overpressure (0.1 MPa) and rock permeability (10 -11 m 2 ) simulated. It is interesting to note that the time delay is more sensitive to rock permeability at low reservoir overpressures. Contrarily, at very high permeability, changes in the gas reservoir overpressure only cause slight changes of the time delay.

Figures 11c andd show the amplitude of the peak in the spectrum for the 12h and 24h components of the simulated CO 2 flux as a function of both rock permeability and gas reservoir overpressure. Some elements of discussion can arise from the comparison between these two figures. In the range of the simulated permeabilities and gas reservoir overpressure, it seems that the CO 2 flux has a 24h component larger than the 12h component, i.e. once again the temperature perturbation seems to be more effective on the simulated degassing. However, there are some regions in the parameters space where the peak amplitude of the 12h component is slightly higher than, or at least comparable with the peak amplitude of the 24h component. This happens for small gas reservoir overpressure (up to 0.01 MPa) and for permeabilities in the range 10 -14 -10 -11 m 2 .

Moreover, the peak amplitudes of the two considered components are comparable for very high permeabilities (∼ 10 -11 m 2 ). In the 12h peak amplitude diagram (Fig. 11c), the iso-amplitude lines are parallel to the overpressure axis, and indicate that this amplitude

D R A F T August 30, 2012, 9:29am D R A F T
is mainly influenced by permeability rather than source pressure changes. A similar, but less marked, pattern is observed in Figure 11d for the amplitude of the 24h component.

Analyzing the changes in 24h and 12h peak in the spectrum, a change in soil permeability will provoke a remarkable and similar increase of both 12h and 24h components, while an increase in the gas reservoir overpressure will induce a different behavior of the two components, with the 24h component increasing while the 12h component decreases (or remains constant). Thus the observation of a correlated change in the 12h and 24h peak amplitudes of the gas flux spectrum may suggest a permeability change in the soil, rather than a change in source conditions as a driving mechanism. It is of note that for low permeabilities, the evolution of the 12h component time delay can also help to distinguish a permeability change from a pressure change.

The described trend could be affected by our choice of initial conditions. The simulations described above were all run with an initial low temperature gradient (3 • C/m).

Additional simulations were performed to asses the role of the initial temperature distribution. Results, shown in Figure 12, suggest that the initial temperature gradient does not affect (or slightly affects) the resulting CO 2 flux. The peak amplitude of both 12h and 24h components is constant in the range 3-30 • C/m (Fig. 12a), and the same is observed for the time delay for both components (Fig. 12b).

Discussion and Conclusion

This work focused on the effects of the atmospheric conditions on soil diffuse degassing. As we are dealing with uprising fluids, in this case higher density hinders the upward motion of the fluid. However, under the conditions considered here, the gas density is very small (on the order of 2 kg/m 3 ) and this term is negligible with respect to the pressure gradient (the pressure difference across the domain is of the order of a few thousands Pa).

The atmospheric pressure acts in different and opposite ways. On one hand, it controls the pressure gradient. In our degassing system, the pore pressure at depth is larger than

D R A F T August 30, 2012, 9:29am D R A F T
at the top. If the atmospheric pressure further diminishes, the gradient across the surface increases, leading to a higher gas flux. Pressure also acts on the CO 2 density, which in this case is directly proportional to it. A high atmospheric pressure act on gas flux reducing the pressure gradient, but increasing the contribution due to density.

The overall value of the gas flux at the surface depends on the complex interplay of all these different contributions, and on their relative magnitude. Under the conditions considered here, the largest changes in the simulated gas flux are associated with the changes of CO 2 mobility, which is mostly controlled by temperature. Pressure effects arise through the pressure gradient term, but affect the degassing to a smaller extent. In this particular case, temperature and pressure are both inversely proportional to gas flux, and their effects can sum up at particular times, when they both follow similar increasing or declining trends.

Applied perturbations and simulated gas flow changes are not only confined at the ground surface, but propagate down to the bottom of the simulated domain. This effect is mainly controlled by the chosen domain permeability and gas reservoir overpressure. The perturbation applied at the surface as described above may be balanced before reaching the reservoir at depth when the fluid flow is slow (as in the case of very low permeability), or affect only the upper part of the domain when the pressure changes applied at the surface are of the same order of magnitude or greater than the overpressure imposed at the reservoir.

Degassing may be affected by many parameters, among which rainfall, air humidity, and wind speed, and then it is not possible to directly compare the absolute values of simulated and observed CO 2 flux. However, results for a base case show a variation in degassing of

D R A F T August 30, 2012, 9:29am D R A F T
the same order of magnitude of the observed changes at Furnas. We calculated two time delays: the correlation between the 24h component of the CO 2 flux with respect the air temperature, and between the 12h component of the gas flux with respect the barometric pressure. Results highlight a delay of ∼ 15 hours for the 24h component, and a delay of ∼ 6 hours for the 12h component, well in agreement with what observed in the real field, although, even in this case, time delays may be affected by other perturbations.

CO 2 flux time series depend on both the domain permeability and the overpressure in the gas reservoir. In order to understand the role played by the rock permeability and by an overpressure driving gas ascent, a dual parameters study was performed. The dual parameters study highlights changes in time delay and peak amplitude for both the considered components of the CO 2 flux (12h and 24h).

Analyzing the time delay, on one hand a mostly constant delay around 6 hours results

for the 12h component, and this delay slightly changes only for very low permeabilities (up to 8 hours); on the other hand the delay for the 24h component presents remarkable changes, ranging from 13 hours for high permeability and overpressure, to 21 hours when both the rock permeability and the gas reservoir overpressure are low.

Amplitude of the peak in the spectrum for both components also changes as a function of rock permeability and gas reservoir overpressure. This study highlights how the CO 2 flux main period depends on the system properties. From the analysis of the peak amplitude, the resulting degassing has mainly a 12h period for low gas reservoir overpressure (up to 0.01 MPa) and for permeability higher than 10 -14 m 2 . The peak amplitude for the two components are comparable for very high permeabilities, and the resulting degassing will present in that case both the 12h and the 24h period.

D R A F T August 30, 2012, 9:29am D R A F T
These results indicate that a coupled analysis of these observables may provide useful hints to discriminate the parameter that is changing within the system. According to our model, the observation of correlated changes in the diurnal and semi-diurnal components of diffuse degassing, under the conditions described above, would be related to a permeability change in the soil. On the contrary, an unrelated evolution of these components would indicate a change in the source conditions. This is in agreement with the observations of [START_REF] Richon | Radon anomaly in the soil of Taal volcano, the Philippines: A likely precursor of the M7.1 Mindoro earthquake[END_REF] that correlated changes in the 12h cycle of Rn time series with a magnitude 7.1 earthquake close to Taal volcano (Philippines). Permeability changes in a volcano/hydrothermal system can be generated by transient stresses, including distant earthquakes [START_REF] Elkhoury | Seismic waves increase permeability[END_REF][START_REF] Manga | Changes in permeability caused by transient stresses: Field observations, experiments, and mechanisms[END_REF], or may result from micro-fracturing associated with intrusive processes and mechanical disturbances, such as a dome collapse, or from mineral dissolution or precipitation [START_REF] White | Permeability Changes During the Evolution of a Geothermal Field Due to the Dissolution and Precipitation of Quartz[END_REF]].

The monitoring of harmonic components of the gas flux signal should thus be integrated in volcano monitoring as it represents a tool to understand the reasons of gas flux changes at the surface.

Acknowledgments. Governing equations solved in TOUGH2/EOS2 for single-phase, singlecomponent (CO 2 ) non-isothermal cases. Symbols: V volume (m 3 ), M mass accumulation term (kg m -3 ), Γ surface area (m 2 ), F Darcy flux vector (kg m 2 s -1 ), n outward unit normal vector, q V volumetric source term (kg m -3 s -1 ), φ porosity, ρ and ρ R fluid and rock density (kg m -3 ), C R heat capacity of the rock formation (J kg -1 K -1 ), T temperature, u CO 2 internal energy (J kg -1 ), κ permeability (m 2 ), µ dynamic viscosity (kg m -1 s -1 ), P total pressure (Pa), λ thermal conductivity (J s -1 m -1 K -1 ), h CO 2 enthalpy (J kg -1 ).

Equations are generally solved for each component and the heat "component". In this case we only consider a component, then the heat is considered as second component

  study: Influence of atmospheric conditions on the CO 2 degassing at Furnas volcano Secondary manifestations of volcanism in the Azores archipelago (Fig. 1) include low temperature fumaroles (maximum temperature around 100 • C), hot and cold CO 2 -rich springs, and several diffuse degassing areas. Continuous monitoring of hydrothermal soil CO 2 flux begun on Furnas volcano, S. Miguel island, in October 2001, with the installation of a permanent gas station that incorporates several meteorological sensors. Daily and seasonal cycles have been observed in the time series of soil CO 2 flux and are coincident

  Two permanent soil CO 2 flux monitoring stations are presently installed inside Furnas volcano caldera (S. Miguel Island, Azores archipelago) (Fig. 1). GFUR1 station, which was running between October 2001 and July 2006, was placed in a garden of the Furnas D R A F T August 30, 2012, 9:29am D R A F T thermal baths, close to Furnas village fumarolic field. Average values measured of CO 2 flux were ∼260 g m -2 d -1 and soil temperature was about 17 • C. The station was reinstalled (and renamed as GFUR3 ) in January 2008 closer to Furnas village main fumaroles, in a thermally anomalous zone where the average soil temperature is 38 • C at 30 cm depth, with a soil CO 2 flux ∼650 g m -2 d -1 . A soil CO 2 flux station (named GFUR2 ) was also installed inside Furnas caldera in the vicinity of Furnas Lake fumarolic ground in October 2004, where soil CO 2 flux values around 350 g m -2 d -1 were measured. Soil temperature at this site is ∼22 • C at about 30 cm of depth. According to Viveiros et al.

Figure 4

 4 Figure 4 describes the 1D domain, 1 m deep, with the rock properties. The domain was

Figure 5

 5 Figure5shows the results of the simulation on the week-long record. The simulated

(

  see Tab. 1). Note that the observed delay at Furnas were obtained by processing CO 2 flux data during several summers and not only during a week, as in our simulation study.D R A F TAugust 30, 2012, 9:29am D R A F T 5. Effect of permeability and source pressure on 12h and 24h CO 2 flux componentsIn this section, a two parameters study will be presented. The changes in peak amplitude in the spectrum and time delay of the 24h and 12h components of the flux variations will be analyzed as a function of the rock permeability and the gas reservoir overpressure. Domain permeability ranges between 10 -15 m 2 and 10 -11 m 2 , and gas reservoir overpressure between 2•10 -5 and 0.1 MPa.

  Harmonic oscillations are observed on the soil CO 2 flux time series at Furnas volcano in S. Miguel island in the Azores archipelago, and these variations are correlated with D R A F T August 30, 2012, 9:29am D R A F T the meteorological variables monitored by permanent stations (wind velocity, barometric pressure and air temperature). Numerical results performed with the TOUGH2/EOS2 simulator show that the CO 2 flux changes when air temperature and barometric pressure change, and presents a time delay with respect to the applied perturbations in agreement with data observed at Furnas. The simulated degassing presents both diurnal and semi-diurnal variations, as effects of applied temperature and pressure, respectively. Lower fluxes are observed when temperature and pressure are higher, in agreement with data from Furnas, where lower degassing is observed during the afternoon. The variable conditions imposed along the upper boundary of the domain affect the gas flux acting on the terms of the Darcy's flow equation: the pressure gradient and the fluid properties. Air temperature acts on fluid mobility (i.e. the ρ/µ term, see Tab. 2): lower temperatures correspond to higher CO 2 density values and lower viscosity. Since the gas flow is directly proportional to its density and inversely proportional to its viscosity (eq. 4 in Tab. 2), a temperature drop has the overall effect of increasing the fluid flow, and viceversa. The gas density also enters in the buoyancy term in the Darcy's flow equation (∇P -ρg).

(dtFigure 2 .Figure 3 .Figure 4 .Figure 5 .Figure 6 .Figure 7 .Figure 8 .Figure 9 .Figure 10 .Figure 11 .Figure 12 .

 23456789101112 Figure 1. (a) Azores archipelago location highlighting S. Miguel Island; (b) Digital

Table 1

 1

	also reports the time

Table 1 .

 1 This work was carried out when the first author was enrolled in the Doctorate school at the University of Bologna, on a grant funded by the Istituto Nazionale di Geofisica e Vulcanologia -Sezione di Bologna. A. P. Rinaldi is currently supported by DOE-LBNL contract number DE-AC02-05CH11231. F. Viveiros is supported by a Postdoc grant from Fundo Regional da Ciência e Tecnologia. Correlation and time delay between each of the meteorological variables and the soil CO 2 flux at GFUR2 for the summer period.

	D R A F T	August 30, 2012, 9:29am	D R A F T

Table 2 .

 2


Frédéric Giroire

Stéphane Pérennes

Issam Tahiri

On the Hardness of Equal Shortest Path Routing

Keywords: OSPF-ECMP, max flow, NP-Hard, approximation OSPF-ECMP, max flot, NP-Hard, approximation

In telecommunication networks packets are carried from a source s to a destination t on a path that is determined by the underlying routing protocol. Most routing protocols belong to the class of shortest path routing protocols. In such protocols, the network operator assigns a length to each link. A packet going from s to t follows a shortest path according to these lengths. For better protection and efficiency, one wishes to use multiple (shortest) paths between two nodes. Therefore the routing protocol must determine how the traffic from s to t is distributed among the shortest paths. In the protocol called OSPF-ECMP (for Open Shortest Path First -Equal Cost Multiple Path) the traffic incoming at every node is uniformly balanced on all outgoing links that are on shortest paths. In that context, the operator task is to determine the "best" link lengths, toward a goal such as maximizing the network throughput for given link capacities.

In this work, we show that the problem of maximizing even a single commodity flow for the OSPF-ECMP protocol cannot be approximated within any constant factor ratio. Besides this main theorem, we derive some positive results which include polynomial-time approximations and an exponential-time exact algorithm. We also prove that despite their weakness, our approximation and exact algorithms are, in a sense, the best possible.

Introduction

In this paper we study the complexity of routing according to the Open Shortest Path First protocol (OSPF) [START_REF]Ospf version 2[END_REF].

The goal of a routing protocol is to make every router capable of deciding, whenever it receives a packet, the next hop router. The decision must be taken locally and quickly, while allowing an efficient usage of the network resources. When OSPF is the routing protocol used in the network, all packets follow shortest paths, according to the link lengths that are set by the network administrator.

When there are several shortest paths between two nodes u and v, the routing depends on the rule that is used for load balancing the traffic among the shortest paths. There are several rules. One of the most used is ECMP (Equal Cost Multiple Path). According to this rule, a router which has several outgoing links on shortest paths toward a destination v, balances the incoming traffic directed to v evenly among all of them.

In order to understand the algorithmic difficulty of OSPF routing, we look into the most essential problem in which the goal is to maximize the throughput when only a single pair of nodes is communicating over the network. We call this problem Max-Ospf-Flow:

• INSTANCE: a directed capacitated graph D = (V, A, c) where V is the set of nodes, A is the set of arcs, and c is the capacity function that assigns to each arc (u, v) a capacity c(u, v), and two vertices s, t ∈ V , respectively the source and the sink of the flow.

• QUESTION: find the maximum flow going from s to t along with its corresponding arc-length assignment under the protocol "Open Shortest path First" using the strategy "Equal Cost Multiple Path".

We give an example of maximum ospf-flow in Figure 1 and a formal description of the problem in Section 2. Let us point out that there may be a large gap between the Max-Ospf-Flow and the standard maximum flow: Figure 2 shows a flow graph with n+2 vertices in which the value of the Max-Ospf-Flow is 2 while the value of the standard maximum flow is n. To the best of our knowledge, this problem has been only proved to be NP-Hard [START_REF] Pióro | On open shortest path first related network optimisation problems[END_REF].

In this paper, we show that the Max-Ospf-Flow cannot be approximated within any constant factor ratio. Nevertheless, we derive several positive results. First, the problem can be ρ-approximated. This implies that, when the capacities of all arcs are equal, the problem can even be solved exactly in polynomial time. We then present an exact exponential algorithm of complexity Õ(2 |A| ) * and show that any exact algorithm is likely to cost 2 Θ(|A|) . Finally, we provide an approximation for networks with small longest distance of factor 1 2 H( cmax cmin ) 1-L † , where L is the length of a longest simple path between the source s and the destination t of the flow.

* Õ(g(n)) = O(g(n) • log k g(n)), for k ∈ N. Logarithmic factors are ignored. † Here H(K) = K i=1 1 i = O(log(K))
denotes the harmonic number of order K.

RR n°8175

Figure 2: A flow graph with a big gap between Max-Ospf-Flow and the standard maximum flow. Arcs on the path s -u 1 -u 2 -• • • -u n have all capacities equal to n while arcs entering t have all capacities equal to 1.

This study shows that the difficulty of Max-Ospf-Flow highly depends on the range of the capacity function and of the "depth" of the network.

Related Work

Some previous works [START_REF] Fortz | Increasing internet capacity using local search[END_REF][START_REF] Izmailov | Administrative weight allocation for pnni routing algorithms[END_REF][START_REF] Bley | On the approximability of the minimum congestion unsplittable shortest path routing problem[END_REF] showed that finding an OSPF routing that optimizes some criterias (load, capacity, . . . ) of the network performance is a difficult task; and this holds for both unsplittable and splittable (with ECMP rule) routings. When the traffic is unsplittable, the lengths have to be set in a way that for every pair of nodes there is a unique shortest path. In [START_REF] Fortz | Increasing internet capacity using local search[END_REF], the authors studied the problem of optimizing OSPF-ECMP lengths so as to minimize the total cost of the network when the cost function on arcs is convex and increasing with the congestion. They showed that this problem is NP-Hard. In particular, they defined the maximum utilization of the network as max a∈A l(a) c(a) where l(a) is the load of arc a, and they proved that minimizing the maximum utilization is NP-Hard. They also provided worstcase results about the performance of OSPF-ECMP routing vs. an optimal multi-commodity flow routing in term of congestion. In [START_REF] Bley | On the approximability of the minimum congestion unsplittable shortest path routing problem[END_REF], the authors studied the unsplittable OSPF that minimizes the maximum utilization of the network. They showed that this problem is hard to approximate within a factor of

O(|V | 1-ε ).
As these problems are difficult, the methods to solve them in practice usually follow a two-phase approach based on linear and integer programming [START_REF] Bley | Integer programming approaches to access and backbone IP-network planning[END_REF][START_REF] Holmberg | Optimization of internet protocol network design and routing[END_REF][START_REF] Dzida | Valid inequalities for a shortest-path routing optimization problem[END_REF]: In the first phase a routing that is not necessarily feasible with OSPF is found and in the second phase lengths that achieve this routing through OSPF are computed, when possible. This second phase attempts to solve what is called the inverse shortest path problem (ISP). This problem can be formulated as a linear program and therefore can be solved in polynomial time. Even if the weights are constrained to take non-negative integer values, ISP remains polynomial thanks to a rounding scheme presented in [START_REF] Ben-Ameur | Internet routing and related topology issues[END_REF]. However finding a solution to ISP that minimizes the maximum length over all arcs is NP-Hard [START_REF] Bley | Inapproximability results for the inverse shortest paths problem with integer lengths and unique shortest paths[END_REF].

In this paper, we consider the problem of approximating the maximum ospf-flow for a network with a single commodity. This problem has been proved to be NP-Hard in [START_REF] Pióro | On open shortest path first related network optimisation problems[END_REF], using a reduction to set cover. To the best of our knowledge, no approximation with non trivial guaranteed ratio has been suggested before for solving this problem.

To conclude, note that, more generally, the results presented here are valid with minor modifications for other shortest path routing protocols as IS-IS [START_REF] Fortz | Optimizing ospf/is-is weights in a changing world[END_REF] or PNNI [START_REF] Izmailov | Administrative weight allocation for pnni routing algorithms[END_REF]. For more information, a survey on the optimization of OSPF routing can be found in [START_REF] Bley | Optimization of OSPF routing in IP networks[END_REF].

Inria 2 Problem Definition

In this section, we show that in the context of our study, which assumes a single commodity, the inverse shortest path problem is guaranteed, under some conditions on the routing, to have a solution and can be easily solved. Therefore, even though Max-Ospf-Flow remains difficult, we can focus on finding an adequate routing function while forgetting about arc-length assignment.

In the case of a single commodity flow, solving the inverse shortest path problem is trivial for any flow f going from s to t. For this we define A = {a ∈ A | f (a) > 0}, and set ∀a ∈ A , l(a) = 0 and ∀a ∈ A , l(a) = 1. Then the shortest paths from the source to the sink are all the paths using only arcs in A . One may wish for a better length function, such that ∀a ∈ A, l(a) > 0, but it exists iff A is acyclic. First, if A contains a cycle then l(a) must be null on all the arcs of the cycle. Second, if A is acyclic, one easily gets a length function since there exists then a potential function p : V → R such that p(u) < p(v) if there is a path from u to v (this potential function is given by the dual problem, see [START_REF] Cook | Combinatorial optimization[END_REF]). Then, one sets ∀a = (u, v) ∈ A, l(u, v) = p(v) -p(u). Notice that all paths from s to t have length p(t) -p(s). Combining those two facts, one can hence always define lengths such that ∀(u, v) ∈ V 2 , all the paths from u to v in A have the same length and with l(a) > 0 on all arcs a not contained in a cycle. Note that, even if it may seem strange, there exist digraphs for which any solution of Max-Ospf-Flow contains a flow loop (see Figure 3). (c) flow with loop Figure 3: Example of an instance where the optimal flow can only be achieved using a loop: when avoiding the loop, one can easily check that the limitation comes from the constraint f /4 ≤ 4 which means that the flow sent from s to t is at most 16; while using a loop the same constraint allows g to be equal to 24 which leads to a flow going from s toward t that is equal to f = g -g/6 = 20. Other configurations with less arcs are not depicted but lead to even lower ospf-flow value.

Therefore, for any flow function, there always exist lengths for which all the paths used by the flow are shortest paths, and this lengths can be non-negative if the flow is acyclic. Nevertheless, it does not mean that any flow function can be achieved using OSPF-ECMP. Indeed f must also fulfill the "Equally Balanced Condition". This means that we must have: at any node a ∈ V and

∀(u, v) ∈ A, either f (u, v) = 0 or f (u, v) = out(u)
where out(u) depends only on u. This allows us to reformulate the Max-Ospf-Flow in a simpler way, without using metric explicitly. But before we need to introduce a few definitions.

Definition 1 (flow graph) Given a directed capacitated graph D = (V, A, c) and two particular nodes RR n°8175 s and t (to be respectively seen as the source and the sink), their corresponding flow graph is defined as the 5-uplet D = (V, A, c, s, t).

Definition 2 (flow function/value) Given a flow graph D = (V, A, c, s, t), a function f : A → R + is a flow function if ∀a ∈ A, f (a) ≤ c(a) and ∀v ∈ V \ {s, t}, (u,v)∈A f (u, v) = (v,u)∈A f (v, u). The flow value is val(f ) = (s,v)∈A f (s, v) -(v,s)∈A f (v, s).
Definition 3 (regular function) A function over a set S is said to be regular if it takes only two values, 0 and another possible one.

Definition 4 (balanced function) Given graph G = (V, A) and V ⊂ V , a function f is said to be balanced on V iff for any u ∈ V the function f u , defined by ∀(u, v) ∈ A, f u (v) = f (u, v), is regular. Definition 5 (ospf-flow function) For a flow graph D = (V, A, c, s, t), a function f : A → R + is an ospf-flow if it is a flow function and if it is balanced on V .
Now, we are ready to give an equivalent formulation of the problem Max-Ospf-Flow:

• Instance : a flow graph D = (V, A, c, s, t).

• Question : find the maximum flow value of an ospf-flow function in A.

The metric problem has only apparently vanished, since choosing on which arc the flow is null and on which arcs it is strictly positive actually determines the underlying hidden length function.

Some notations For f a function from a set S to R, we denote f (S) = s∈S f (s). Given a digraph with vertex set V and arc set A we also abusively denote (U, V ) the set of couples {(u, v) ∈ A | u ∈ U, v ∈ V }. We also note the set of integers {1, . . . , N } by [N ].

Inapproximability results

We use a reduction to the MAX-3-SAT problem in order to prove first that Max-Ospf-Flow is hard to approximate, in a polynomial time, within an approximation ratio of 31 32 . Then we use a self-amplifying method to prove that this problem is not in APX. We recall that MAX-3-SAT is the problem defined by:

• INSTANCE : a set F = {B j } j∈[m] of clauses, built on a finite set {x i } i∈[n] of binary variables,
such that each clause contains exactly three literals.

• QUESTION : find a truth assignment of the variables {x i } i∈[n] that maximizes the number of satisfied clauses in F .

Unless P = NP, Problem MAX-3-SAT cannot be approximated within a factor 7 8 + ε for any constant ε > 0, see [START_REF] Håstad | Some optimal inapproximability results[END_REF]. In particular, even when the optimum value is m, it is still NP-Hard to find an assignment that satisfies more than 7 8 + ε clauses of F , ∀ε > 0.

Proposition 1 It is NP-hard to approximate Max-Ospf-Flow within a factor 1 -l-3 8l + ε, ∀ε > 0, even if ∀a ∈ A, c(a) ∈ {1, l}. As example, it is NP-hard to approximate Max-Ospf-Flow within a factor 31 32 + ε, ∀ε > 0, even if ∀a ∈ A, c(a) ∈ {1, 4}, Inria Proof. For seek of simplicity, our reduction uses a multigraph. To get a graph from it one can simply insert a node on every arc. This will make the graph simple and will only double its initial size. Given an instance I of MAX-3-SAT, defined by the set of binary variables {x i } i∈[n] and the set of clauses {B j } j∈[m] , we build the following flow graph D I = (V, A, c, s, t):

-For each i ∈ [n], we create a vertex a i representing variable x i .

-For each j ∈ [m], we create a vertex b j representing clause B j .

-For each i ∈ [n], we add k i parallel arcs with capacity 1 going from s to a i , where k i is the number of times the variable x i appears in the clauses.

-We add an arc with capacity 1 from every b j to t, with j ∈ [m].

In addition to those arcs, we add arcs that depend on the structure of the clauses :

-

∀(i, j) ∈ [n] × [m]
, if B j contains the positive literal x i we add an arc with capacity

1 from a i to b j ; -∀(i, j) ∈ [n] × [m]
, if B j contains a negative literal xi we add l parallel arcs with capacity 1 l from a i to b j An example of the construction is displayed in Figure 4. Note that the size of D I is linear in the number of clauses.

B 1 = (x 1 ∨ x2 ∨ x 4 ), B 2 = (x 2 ∨ x3 ∨ x 4 ) and B 3 = (x 2 ∨ x 3 ∨ x4 ).
[The example assumes l = 4. Dashed arcs have capacity 1 4 .

The others have capacity 1].

The idea of the reduction is to emulate binary variables as follows: For (i, j) ∈ [n] × [m] we define the virtual link vl(a i , b j ) as the set of arcs connecting a i to b j , we say that the virtual link vl(a i , b j ) is positive (resp. negative) virtual link if x i (resp. x i ) is a literal of B j . Then, one must choose at a i either to send significative flow on the positive virtual links or the negative ones, and this is exclusive. This formalizes with the following property:

(a i , b k ) ∀b k , k ∈ [m] is at most 1 l .
To prove (a) remark that positive virtual links are made of a single arc with capacity 1, so if (a) premise happen we must have f (a i , b j ) > 1 l and since f ai is regular no flow can leave a i on the arcs with capacity 1 l . Case (b) is similar, since (b) premise means that the arcs with capacity 1 l are used, and since f ai is regular, so f (a i , b j ) ≤ 1 l (since there is a unique arc for the positive virtual links). We first prove that the Max-Ospf-Flow on D I is m if I has a truth assignment that satisfies all clauses. To do so, we associate to this truth assignment a flow function f on D I as follows: for each satisfied clause B j , we choose one of its satisfied literals, and if it is associated to the variable x i we send one unit of flow on the path s -a i -b j -t using the virtual virtual link corresponding to the literal. The flow f has then a value m and is balanced since:

-∀i ∈ [n] with x i true all arcs leaving a i have flow 0 or 1.

-∀i ∈ [n] with x i false all arcs leaving a i have flow 0 or 1 l . Therefore f is an ospf-flow. Moreover, this flow is maximum since the cut ({t}, V \ {t}) has a capacity m. Consequently, the value of the Max-Ospf-Flow on D I is m.

We now assume that we are able to find (using a polynomial algorithm) an ospf-flow function f . To that ospf-flow flow we associate a binary assignment for each i ∈ [n] as follows: a) If a vertex a i sends strictly more than 1 l units of flow along a positive virtual link we set x i to true. b) If a vertex a i sends strictly more than 1 l units of flow along negative virtual link we set x i to false. According to Property 1 this assignment is consistent. Note that the above rules may let some variables unassigned, in which case we simply discard them (or let them unassigned).

Consider now one of the nodes b j , j ∈ [m] for which the flow f (b i , t) is greater than 3 l . Since there are at most three virtual links entering the node b j one virtual link (a i , b j ), i ∈ [n] carries a flow strictly greater than 1 l . Then, either vl(a i , b j ) is a positive virtual link and x i was set to true (by rule (a)) and the clause B j is satisfied. Or vl(a i , b j ) is a negative virtual link and x i is set to false (by rule (b)) and the clause B j is satisfied. So such a node b j corresponds a clause B j satisfied by our assignment. So (unless P = NP) we have at most ( 7 8 + ε)m such nodes. But then, the flow is at most

7 8 + ε m • 1 + 1 -( 7 8 + ε) m • 3 l = m -l-3 l 1 8 -ε m = 1 -l-3 8l m + ε l-3 l
m. So we conclude that ∀l ≥ 3 and ∀ε , (unless P = NP) no polynomial algorithm computes a flow with value greater than (1 -l-3 8l + ε )m. Note that the constructed graph uses capacities { 1 l , 1} so we multiply all the capacities by l to get an equivalent graph with capacities in {1, l}.

Theorem 1 It is NP-Hard to approximate Max-Ospf-Flow within any constant factor.

Proof. To prove the result, we use a self-amplifying error technic. Indeed, the flow graph D I associated to a 3-SAT instance I connects s to t is like a virtual arc (s, t) with capacity m, but for which it is NP-Hard to use a capacity larger than 31m 32 . Considering a flow graph D, we will replace each of its edges by those virtual arcs. Two factors of error will then get multiplied: one because no polynomial algorithm can use the full capacity of the virtual arcs, the other because there is another error on the original network itself. We note OPT(D) the value of Max-Ospf-Flow on D. Let us introduce two definitions: Inria • For a flow graph D = (V, A, c, s, t), we define x ⊗ D = (V, A, c x , s, t) where c x is defined by ∀a ∈ A, c x (a) = x • c(a). For any positive number x, there is a trivial bijection which relates every flow function f for D to a flow function g for x ⊗ D, defined by ∀e ∈ A, g(a) = x • f (a).

• Given two flow graphs D 1 and D 2 , we define D 1 • D 2 as the graph in which every arc a = (u, v) ∈ A (A is in D 1 ) with capacity c(a) is substituted with the following sequence (see Figure 5): a copy of c(a) ⊗ D 2 with source s(u, v) and sink t(u, v), two arcs (u, s(u, v)) and (t(u, v), v) with capacities c(a)M , where M is an upper bound of OPT(D 2 ). Suppose that one had proven that the Max-Ospf-Flow is hard to approximate within ρ 1 and ρ 2 ; namely we can associate to any instance I of an NP-Complete problem a flow digraph D 1 (resp. D 2 ) such that a ρ 1 (resp. ρ 2 ) approximation of the Max-Ospf-Flow allows to solve I. We build the flow graph D 1 • D 2 . Note first that the size of D 1 • D 2 is the product of D 1 and D 2 sizes (hence if D 1 and D 2 have polynomial size then so does D 1 • D 2 ). The maximum ospf-flow of D 1 • D 2 has value OPT(D 1 )OPT(D 2 ). Suppose now that we can find in polynomial time an ospf-flow of value strictly greater than ρ 1 ρ 2 OPT(D 1 )OPT(D 2 ) for every instance I. For any instance I, two different cases can happen:

• Either, there exists at least one virtual arc with a flow value strictly larger than ρ 2 OPT(D 2 ). In this case, we use reduction 2 and solve the NP-problem for instance I.

• Or, for all virtual arcs, the flow is less than ρ 2 OPT(D 2 ). In this case, we can define a valid ospfflow for the graph OPT(D 2 ) ⊗ D 1 by simply taking as flow value of arc (u, v) the flow value of arc

(u, s(u, v)) in D 1 • D 2 .
We can build a valid ospf-flow for D 1 by dividing the value of the flow on each arc by OPT(D 2 ). This ospf-flow is of value greater than ρ 1 OPT(D 1 ). We then use reduction 1 and we solve the NP-problem for instance I.

Hence, by contradiction, it is not possible to find in polynomial time an ospf-flow of value greater than ρ 1 ρ 2 OPT(D 1 )OPT(D 2 ) for every instance I. To complete the proof we start with ρ 1 = ρ 2 = ρ = 31 32 . And starting from hardness to approximate within ρ we prove hardness to approximate within ρ 2 it follows that Max-Ospf-Flow cannot be approximated within any constant unless P = NP.

Positive results

Even though it is hard to find the optimal solution to Max-Ospf-Flow, we can still get an approximation algorithm.

Theorem 2 Given a flow graph D = (V, A, c, s, t). Let c min and c max be respectively the smallest and the largest value of the capacity function c. Then Max-Ospf-Flow can be approximated on D in polynomial time within a factor c min /c max .

Proof. Let OPT be the maximum value of an ospf-flow on D, and let D i to be the directed graph that have the same elements as D but where the capacity function is the constant function that assigns i to every arc. We compute a maximum integral flow function F 1 on D 1 . On the one hand, the value of an ospf-flow on D cannot exceed the maximum value of a flow on D cmax , which is c max v(F 1 ), so OPT ≤ c max v(F 1 ). On the other hand, c min F 1 is a feasible ospf-flow for D with value c min v(F 1 ) ≥ cmin cmax OPT. When c min = c max we get the next corollary :

Corollary 1 Any instance of Max-Ospf-Flow where all capacities of the flow graph are equal can be solved in polynomial-time.

We now give an exact algorithm with exponentional solving time. Proof. Consider an ospf-flow function f x with a flow value x and remark that the knowledge of the configuration A = {a ∈ A | f x (a) > 0} entirely determines f x in the following sense: For a given configuration, there exists at most one ospf-flow of value x. Indeed if d(v) denotes the out-degree in A of a vertex v the following equations are satisfied:

∀v ∈ {s, t}, ∀(v, w) ∈ A f x (v, w) = 1 d(v) (u,v)∈A f x (u, v) ∀(s, w) ∈ A f x (s, w) = x d(s)
Note that the above system is well-defined iff all nodes adjacent to some edge in A belong to a path in A from s to t. Since the configuration associated to the Max-Ospf-Flow satisfies this condition, we can discard any degenerated set A violating the above condition.

If f 1 is the function that solves the above system for x = 1, then for x ∈ R the solution is given by the function f x = x • f 1 . Therefore, the maximum flow value of an ospf-flow using the configuration A , noted val(A ), expresses as the maximum among all x ∈ R that respects the following capacity constraints:

∀(v, w) ∈ A , xf 1 (v, w) ≤ c(v, w)
Hence, we have val(A ) = min (v,w)∈A c(v,w) f1(v,w) . Finally, to compute the max-ospf-flow, we can simply compute val(A ), in polynomial time, for each non degenerated configuration and take the maximum.

Remark 1 In the above proof, when A forms an acyclic graph, the system is diagonal and can be solved very efficiently by pushing the flow values from the source node. But, as mentionned in the introduction, we cannot limit our search to such acyclic configurations since there exist flow graphs for which the maximum ospf-flow contains cycles.

We now give a more efficient approximation algorithm when the range of capacities is large. The algorithm is based on a sampling argument given in Lemma 1 and on the monotonicity of a (multisource) ospf-flow, meaning that if a ospf-flow exists, ospf-flows of smaller values always exist in the sense of Lemma 2.

Definition 6 Given two functions f, g from S → R, we say that g is f -dominated if ∀x ∈ S, f (x) ≤ g(x). ‡ Õ(g(n)) = O(g(n) • log k g(n))

Inria

Lemma 1 Let S be a finite set, and a function p : S → [M AX], then there exists an integral regular p-dominated function q with q(S) ≥ p(S) H(M AX) § .

Proof. Assume that the unique non zero value that we pick for q is i. Then, since q ≤ p, the best choice for q (in order to maximize q(S)) is to simply set q(x) = 0 iff p(x) < i and q(x) = i iff p(x) ≥ i. So, let us define the two sets: ∆(i) = p -1 (i) and G(i) = ∪ j∈[M AX],j≥i ∆(j). Then, from definition, for any i ∈ [M AX], we can build a regular function q, p-dominated, such that q(S) = i • |G(i)|. It follows that the best value that we can get for q(S) is:

OPT = max i∈[M AX] i • |G(i)|.
Note also that, with those notations (using a standard double couting argument), we have:

p(S) = i∈[M AX] i|∆(i)| = i∈[M AX] |G(i)|. We have ∀i ∈ [M AX], |G(i)| ≤ OPT i . So p(S) = AX] |G(i)| ≤ OPT i∈[M AX] 1 i = H(M AX)OPT. It follows that OPT ≥ p(S)
H(M AX) . Note that by construction the optimal function q is integral and OPT is integral, so we indeed have OPT ≥ p(S)

H(M AX) .
Definition 7 (Multisource ospf-flow) Consider a digraph G = (V, A) and a capacity function on the arcs c, a multisource ospf-flow with n sources {s i } i∈[n] and a sink t is a function f :

A → R + such that • ∀a ∈ A, f (a) ≤ c(a) and • ∀v ∈ V \ {s i } i∈[n] ∪ {t}, uv∈A f (u, v) = vu∈A f (v, u). • balanced on V . For i ∈ [n], val i = (si,v)∈A f (s i , v) -(v,si)∈A f (v, s i ) is called the ouput value of source i The flow value of f is defined as val(f ) = i∈[n] val i .
Lemma 2 (monotonicity) Consider a flow graph with k sources {s i } i∈[k] and a sink t. If it exists a multisource ospf-flow in which the source s i has an output value val i , then, for any {val i } i∈[k] such that ∀i ∈ [k], val i ≤ val i , there exists a multisource ospf-flow in which each source has an output value val i .

Proof. We take the configuration A associated to the initial ospf-flow, i.e. A = {a ∈ A; f (a) > 0} . We then want to determine if an ospf-flow with output values val i , i ∈ [k] on the sources is feasible with configuration A . Such a flow is entirely determined by the values val i , i ∈ [k] (see the proof of Theorem 3), and the flow on an arc a is of the form i∈[k] λ a,i val i where λ a,i are all positive or null. So the load of an arc is increasing with the source outputs and the result follows (the set of feasible ospfflows, with configuration A , is defined by ∀a ∈ A, λ a,i val i ≤ c(a)). Definition 8 For l ∈ N, we say that a flow graph D = (V, A, c, s, t) is l-layered if its vertex set can be partitionned into

V = ∪ i∈[l] val i such that V 1 = {s}, V l = {t} and the arc set A satisfies A = ∪ i∈[l-1] A i where ∀i ∈ [l -1], A i ⊆ (val i , V i+1 ).
Proposition 2 Given an integer l ≥ 3, let D = (V, A, c, s, t) be an l-layered flow graph. Given an integral flow function F that has values in [c], we can compute in polynomial-time an ospf-flow of value val(f ) greater than H(c) -l val(F ).

Proof. The idea is to sample the original flow function F level per level. Let ρ = H(c) -1 be our sampling factor. We build by induction on t a flow function g t that has the following properties:

g t is balanced on the t last levels (i.e. on the vertex set

B t = ∪ i∈[l-t,l] V i ). -g t is F -dominated, that is ∀a ∈ A, g t (a) ≤ F (a). -v(g t ) ≥ ρ t • val(F ).
g t is integral till level t (i.e ∀a ∈ ∪ i∈[t] A i , g t (a) is an integer).

For t = 0, we take g 0 = F . Since B 0 is reduced to the sink, the inductive hypothesis holds.

We now build g t from g t-1 for t ≥ 1. We start setting g t = g t-1 . By induction hypothesis, g t is balanced on all the levels in [l -t + 1, l] (i.e on B t-1 ). We want to make g t balanced on level V l-t . So, we look at the situation of the outgoing flow from level V l-t to level V l-t+1 . We consider the flow function g t-1 on the set A l-t . By induction hypothesis, g t-1 is integral with values in [c]. So we can sample g t-1 on the set A l-t using Lemma 1 and get an integral, regular, g t-1 dominated function q t such that q(A l-t ) ≥ ρg t-1 (A l-t ). We then set ∀a ∈ A l-t , g t (a) = q t (a). Hence, the flow from level l -t to level l -t + 1 is at least ρv(g t-1 ) = ρ t val(F ).

The function g t is not anymore a flow since it does not respect the conservation anymore on levels V l-t and V l-t+1 . Nodes in V l-t receive too much flow, and nodes in V l-t+1 send too much flow.

-For nodes in V l-t+1 , we use monotony (Lemma 2) which allows us to decrease the flow sent toward the sink while keeping the ospf constraints satisfied.

-For nodes in V l-t we can easily decrease the flow since the flow function is not required to fulfill the regularity condition. Nevertheless, we have to keep the flow integral. But since q t is an integral function, the excess flow received at any vertex v ∈ V t-l is an integer. So we can decrease the flow unit per unit and keep g t integral on all levels till t -l.

So all the inductive conditions are fulfilled for t, the result follows by induction.

We define L D (u, v) as the length of the longest simple path between u and v in a flow graph D.

Theorem 4 Let D = (V, A, c, s, t) be a flow graph, and assume that c : A → [c], then Max-Ospf-Flow can be approximated in polynomial time within a factor H(c) 1-L D (s,t)

Proof. We first compute a maximum integral flow F without flow loops (such always exists). The arcs used by the flow then form an acyclic digraph, and we only consider this digraph. We now assign to each vertex a level as follows: at step i ≥ 0, we mark any vertex for which all the in-neighbors were marked at previous steps and assign to such a vertex the Level i. The procedure ensures that the arcs all go from a level to a greater level. Note that the source is the only vertex at Level 0 and that the sink is the only vertex at Level L D (s, t), since there is a path going from the sink to all vertices and that all paths go to the sink. So we have almost layered our digraph, to complete the process we add k -1 intermediary nodes on any arc going from level i to level i + k in order to get only arcs between consecutive levels. See Figure 6 for an example.

Then we apply Proposition 2 to the layered digraph with l = L D (s, t) + 1. We obtain in polynomial time an ospf-flow f on the layered digraph of value greater than H(c) -l val(F ). We can easily build an Inria ospf-flow of same value on D (the flow on all arcs of a path with intermediary nodes is equal, we assigne this flow value to the corresponding arc of the original digraph). This ends the proof since the value of a maximum flow F is an upper bound on the value of a maximum ospf-flow.

a s t F 4 F 5 F 3 F 2 F 1 a b s t F 3 F 2 F 1 F 4 F 2 F 1 F 5 F 1

Conclusion & Open problems

In this study, our main incentive was to investigate the difficulty of routing the data "optimally" in networks that have opted for OSPF with ECMP strategy, as a dynamic routing protocol. To do so, we concentrated on the fondamental problem of maximizing the flow from a single source to a single destination on a network where this protocol is assumed to be running. The main result shows that there exists no constant factor approximation if the capacities are not bounded. When all capacities are integers in [c], There is a c-factor approximation, however, we have been unable to find an algorithm with ratio better than c So we conjecture that for any µ Max-Ospf-Flow cannot be approximated within a factor better than c 1-µ . Nevertheless, if the longest path from s to t is bounded by L, then we derive an approximation algorithm with the ratio 1 2 H( cmax cmin ) 1-L . Additionally, we suggest an exact algorithm that has a complexity of Õ(2 |A| ). It is an open problem to find an algorithm with a better exponential base. However, there is a limitation on the efficiency of such algorithm, since it would be able to solve MAX-3-SAT (see the reduction of Theorem 1) and the best known algorithm for solving an instance with m clauses has complexity Õ(1.3 m ).

For the proof of Theorem 5, we use the notion of multi-source ospf-flow. In fact, most of the results of this paper could be extended to this case, in which there still is a single commodity, but in which the flow can be sent by several sources and, potentially, also sent to several destinations. This case could

Figure 1 :

 1 Figure 1: Example of maximum ospf-flow of value 14.

Figure 4 :

 4 Figure 4: Flow graph D I constructed from a set F containing the clausesB 1 = (x 1 ∨ x2 ∨ x 4 ), B 2 = (x 2 ∨ x3 ∨ x 4 ) and B 3 = (x 2 ∨ x 3 ∨ x4 ).[The example assumes l = 4. Dashed arcs have capacity1 4 . The others have capacity 1].

Property 1

 1 Any ospf-flow function must satisfy RR n°8175 a) If for some (i, j) ∈ [n] × [m] the flow on a positive virtual link vl(a i , b j ) is strictly greater than 1 l then any flow on a negative virtual link vl(a i , b k ) ∀b k , k ∈ [m] is at most 0. b) If for some (i, j) ∈ [n] × [m] the flow on a negative virtual link vl(a i , b j ) is strictly positive then any flow on a positive virtual link vl

Figure 5 :

 5 Figure 5: A virtual arc of D 1 • D 2 (Right) substituting the single arc (u, v) of D 1 (Left).

Theorem 3

 3 Given a flow graph D = (V, A, c, s, t), we have an exact algorithm for solving Max-Ospf-Flow on D which runs in Õ(2 |A| ) ‡ .

  b

Figure 6 :Theorem 5 2 cmax cmin f 0 . 2 cmax cmin f 0 .

 652020 Figure 6: Example of the construction of Theorem 5 Left: D . Right: D. A vertex u is placed at the level L D (s, v), with L D (s, v) the length of a longest path between s and v. Each arc (u, v) ∈ A is then subdivided into a path of length L D (u, v).

  2004 route des Lucioles -BP 93 06902 Sophia Antipolis Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr

RR n°8175

This work has been founded by ANR-DIMAGREEN, ANR-AGAPE, ANR-GRATEL and RAISOM

be interesting to study further as it models Content Delivery Networks (CDNs), estimated to account nowadays for 30 to 50% of Internet global traffic. In this context, data is replicated in several data centers and users in different geographical places get the data from these different sources.