Wavelet-based analysis of ground deformation coupling data from satellites (Sentinel-1, SMOS) and from shallow and deep wells in Southwestern France

Acquisitions of the Sentinel-1 satellite are processed and comprehensively analyzed to investigate the ground displacement during a three-year period above a gas storage site in Southwestern France. Despite quite low vertical displacements (between 4 and 8 mm) compared to the noise level, the local displacements reflects the variations due to charge and discharge during summer and winter periods, respectively. A simplified mechanical model can explain these displacements at both storage sites (Lussagnet and Izaute) by the gas exploitation. However, ground movements of low-magnitude may be induced by natural factors, such as the temperature or the surface soil moisture (SSM). Using an additive decomposition, we show first that the temperature has a second order contribution compared to the reservoir pressure. Using a wavelet-based analysis, we show there is an uplift in the Lussagnet zone that contrasts both in phase and period with the seasonal deformation and that is linked to the SSM measured by the SMOS satellite. This other displacement is consistent with the water infiltration in the unsaturated zone followed by the swelling of a clay layer. This work provides therefore new insights on the ground deformation using a three-year integrated monitoring of a gas storage site.


Introduction
The Underground Gas Storages (UGS) are designed to address different needs that include a strategic gas reserve, the regulation of the gas supply, meeting seasonal peak heating and electricity demand. Gas is stored from spring to autumn when the demand is lower and withdrawn from October to April, when the demand is higher during the winter period. In France, there is moreover an increasing concern to balance daily load needs related to the increase of Subsidence in Mexico city center 3,4 , uplift in Brussels city center 5 and both phenomena in Phoenix 6 are just a few recent of many extensively studies using this remote sensing technique. Conversely, the monitoring of gas storage using the DInSAR techniques has received until now few attentions in scientific literature except for CO2 geological storage 7 . A major limiting factor to this purpose was the non-availability of both spatially and temporally high-resolution Synthetic Aperture Radar (SAR) dataset. Another limiting factor is the low amplitude of deformation associated to gas storage in deep geological layers; the displacement to be measured is in this case of the same order of DInSAR measurement precision (typically a few millimeters). Finally, it has been shown that DInSAR techniques are more performant on regular motions than on displacements that have nonlinear behavior respect time (i.e with strongly varying rates) as such irregular displacements require finer temporal sampling to be better characterized 8  April 2016) that work in a pre-programmed operation mode to produce a consistent long-term data archive built for applications based on long time series 9 . SAR images are acquired in interferometric mode every 6 days (combining Sentinel-1A/1B) over Europe since October 2016.
Due to the potential substantial damage to buildings and infrastructure, the mapping of expansive soils and the quantification of the clay swelling potential is a major concern for natural risk prevention plans. In France, the shrink/swell risk is the second most important cause of financial compensation from insurance companies behind the flood risk. In 2010, based on the spatial distribution of infrastructure damage and the stratigraphy at the resolution of sedimentary basins, the French Geological Survey (BRGM) published a predictive 1:50 000 swelling-risk map of France. This map indexed the territory as (i) no, (ii) low, (iii) moderate, or (iv) high risk. In the studied zone in the Aquitaine Basin, there is a low or a moderate swelling risk. At this resolution, the heterogeneity of the mineralogical composition of the sedimentary formations is not considered [10][11][12] 13 . Another way to improve this swelling-risk map is to use remote sensing satellite or aerial photography. Until now, the use of DInSAR was not operational for the mapping of swelling clays, mainly because of the non-availability of both spatially and temporally high-resolution and high-quality SAR dataset suitable to the very high variability of such surface deformation phenomena. Twenty years later, the question "Can we map swelling clays with remote sensing?" asked by Van der Meer in 1999 14 is still a topical issue.
This work investigates therefore three main questions: (1) Are the variations detected in DInSAR measurements strongly related to the underground gas storage operations at Lussagnet and Izaute sites? ; (2) Can the monitoring based on DInSAR processing be used to monitor the movement of the Gas-Water contact in the reservoir? ; (3) Can the DInSAR processing be used to assess locally the shrink/swell hazard as given by the 1:50 000 geological map? The results section of this paper is organized around these three questions.

Geographical and Geological Setting
The underground gas storage in the studied area is a double storage site lying at the boundaries between both departments "Les Landes" and "Le Gers", situated in the Aquitaine basin in Southwestern France, about 100 km north of the Pyrenees mountain (Fig. 1). The total reservoir structure is an anticline with two culminations, which are only approximately 10 km apart, the Lussagnet reservoir at the west side and the Izaute reservoir at the east side (Fig. 2). The highest point of the top of Lussagnet reservoir is located at a depth of about 545 m below ground level (mbgl) with a thickness of 40 m; the top of Izaute reservoir is at approximatively at 500 mbgl 15 . Both reservoirs are deposited during the Eocene, and are composed of fair consolidated sandstones, called "infra-mollassic sand", with some interlayered claystones 16 . The hydrodynamic parameters of these reservoirs are variable. The total mean porosity of the sandstones varies from 20% to 35%, and their average permeability from 1 to 10 Darcy. The Lussagnet and Izaute gas storage total capacity are respectively equal to 2.9 and 3.0 10 12 m 3 (referred to normal conditions) and only the working gas, respectively 1.4 and 1.5 10 12 m 3 , is stored and withdrawn during the exploitation cycle. The remaining part, called cushion gas, supplies pressure support and prevents surface installations from excessive water production.
Concerning the surface layer, the shrink/swell hazard with 4 levels (null, low, moderate and high) has already been evaluated in both departments "Les Landes" and "Le Gers" using different criteria, including the surface geological map at the scale 1:50 000 and other geotechnical and mineralogical criteria 17,18 . In particular, the Tortonien clay surface layer (called "m5" in the 1:50 000 geological map) is characterized by a moderate shrink/swell hazard and the Serravallien sand surface layer (called "m4") by a low hazard. Just around the gas exploitation wells, there is a moderate shrink/swell hazard at the Lussagnet site and a low hazard at the Izaute site ( Fig. 2 and 4).

Sentinel-1 DInSAR processing
SqueeSAR™ is the proprietary multi-interferogram technique which provides high precision measurements of ground displacement by processing multi-temporal satellite radar images acquired over the same area 19 . The methodology exploits the so called Persistent Scatterers (PS) but also lowamplitude but coherent returns that are identified on a pixel-by-pixel basis called Two major biases can affect SAR interferograms. Firstly, non-stationary tropospheric SAR signal delay can add an undesired DInSAR signal component.
This delay can partly -and when the local topography is significant -be correlated to topography, in which case it can be modeled and removed by linear regression with a digital elevation model 20 . Data stacking is also an effective way to reduce the non-topography-related atmospheric noise. Secondly, random temporal changes on the surface of the Earth can reduce the signal to noise ratio (snr), which is characterized through the Interferometric coherence.  bottomhole pressure, rainfall, surface soil moisture, piezometric level). Usually, they present a shorter time sampling (typically 1 day) and consequently, they must be down-sampled to the revisiting time period of Sentinel-1 (12 days), e.g. for the bottomhole pressure or the piezometric level. However, some other variables require more complex processing than a simple down-sampling. For the rainfall, a cumulated rainfall during the Sentinel-1A period (12 days) is recalculated. For the surface soil moisture, an average value for each 12-day period has been calculated using the 10-day SMOS Level 3 product as explained in the next section.

SMOS Level 3 product (Centre Aval de Traitement des Données)
In situ measurements of soil moisture are limited in their spatial and temporal extent. Satellites provide more extensive spatial coverage and have a temporal resolution ranging from 1 to 35 day(s). We use here the term Surface processing. In order to prevent any inconsistency resulting from interpolation over highly heterogeneous surfaces, no spatial averaging is operated in the algorithms.
It must also be noted that ascending and descending overpasses are bound to show different values of the retrieved parameters that may not be always comparable and they are thus retrieved separately. The SMOS Level 3 SSM (RE04v300) products were downloaded through the website of the Centre Aval de Traitement des Données SMOS (CATDS, https://www.catds.fr/). The data are presented over the Equal-Area Scalable Earth (EASE grid 2) 24 with a sampling of about 25 km x 25 km (Fig. 1). The performance of each satellite SSM product depends on many factors such as, but not limited to, soil type, climate, presence of noise (Radio Frequency Interference) and land cover. It is therefore difficult to predict the performance of SSM products over a region, without performing a quality assessment using in situ measurements on that region.

SqueeSAR TM time series analysis around gas exploitations wells
To investigate the DInSAR products in both gas exploitation zones, eight and fourteen DInSAR measurement points are respectively selected at Lussagnet site (4 PS and 4 DS) and Izaute site (11 PS and 3 DS) (Fig. 4). This selection is mainly based on the measurement points quality parameters. As explained in the

Vertical displacement due to reservoir pressure variations
Qualitatively, from the visual comparison between the mean LOS displacement at both sites and the pressure of both storage reservoirs (Fig. 6), where cylindrical coordinate ( , ) is taken (Fig. 8A), and 0 ( ) is the Bessel function of the first kind for integer order 0. For taking into account of realistic reservoir geometry, we subdivide it by many point sources in order to conserve the reservoir volume and we take a summation over all the contribution of each part of the reservoir. The assumed model parameters are summarized in Table 1.
This is a phenomenological, namely kinematic, expression without solving hydrogeological or geomechanical system of the reservoir. From the geological map (Fig. 2), we estimate the extension of the reservoir as 5.4 km 2 for Lussagnet and 13.1 km 2 for Izaute, respectively (Fig. 8A). The gas is stored at the top of the reservoir such that the aquifer level is lowered. Schematically we simplify this by a vertical cylindrical reservoir in which the injected gas pushes down the native water. The modeling geometry is therefore simply vertical so that the deformation pattern calculated through Equation (3) concentrates around the reservoir surface ( Fig. 8B). From this estimation, the vertical displacement on the ground surface arrives at 8 mm and 4 mm for Lussagnet and Izaute, respectively (Fig. 8C). The deformation pattern is quite independent, meaning that the deformation amplitude is close to zero half way in-between the two reservoirs. The amplitude is consistent with the DInSAR measurements in the Lussagnet case (7.4 ± 1.7 mm), but lower than the DInSAR measurements in the Izaute case (6.4 ± 0.8 mm) (see previous section). As it is clear from Equation (1), the amplitude of deformation depends strongly on the material constants ( , , ) given in Table 1. In the case where these three parameters vary with ±10%, changes between +41% and -30%.
In the Izaute case, this gives an upper bound of the simulated vertical displacement of 5.6 mm, which is consistent with the DInSAR measurement (6.4 ± 0.8 mm). It can also be noted that a smaller reservoir geometry shall lead to a  (Fig. 6A). We do not observe other significant differences between the LOS displacement and the reservoir pressure for other sub-periods in the Lussagnet case and in the Izaute case (Fig. 6B). Our hypothesis is that this significant difference during this short period might be due to the clay swelling in the Lussagnet zone, as supported by the current swelling risk map ( Fig. 2 and 4A,B).
According to our simplified model, the monitoring wells in both sites are outside the domain influenced by the gas injection, i.e. the vertical displacement amplitude due to the gas injection at these monitoring wells is less than 1.25 mm (corresponding to a maximum peak-to-peak less than 2.5 mm) (Fig. 8B).
Therefore, it is sound to investigate the LOS displacement at the monitoring wells around both gas exploitation zones, in order to decouple the vertical displacement due to the gas injection from other potential phenomena (e.g. clay swelling). First, we focus our attention on the comparison between LOS displacements at the aquifer monitoring wells around the gas bubble and piezometric levels during the first 2 years and half ( Fig. 9 and 10). In contrast to the result in the gas exploitation zone itself, there is no obvious correlation between these two time series and the signal-to-noise ratio of the LOS time series is low. To investigate further these LOS displacements, we use wavelet-based tools which have demonstrated their usefulness in this context of seasonal variations to extract features from low signal-to-noise ratio time-series (see the method section). The proposed methodology is to study the time-frequency relationships between the surface displacements measured by DInSAR and three potential natural triggering factors: (1) the cumulated rainfall acquired by the Mont-de-Marsan Weather Station (Fig. 1), (2) the surface soil moisture acquired by the SMOS satellite in a 25km grid cell around the studied area ( Fig. 1), (3) the shallow phreatic groundwater piezometric level measured at Latrille (Fig. 2).

Wavelet-based analysis of DInSAR time series around
Where t is the time in days (relative to 18 October 2014), Amax the amplitude of 2 ± 0.32 mm, ϕ the phase delay of 32 days days and T the period of 365 days. We add to this sinusoidal model the LOS displacement at L2, supposing that the vertical displacement due to the clay swelling is the same at both locations because of the same depth of the clay layer: Indeed, this new model M(t) fits quite well the displacement for the DS called BKB8Z5Y (near L1) between December 2015 and April 2016 (Fig. 17).
In conclusion, the combination of a LOS displacement due to the pressure change and a displacement due to the swelling of a 3-m depth clay layer fits well during the whole monitoring three-year period (Fig. 18C). Additionally, some power signals during the first semester of year 2 (between December 2015 and April 2016) with a significant level against red noise can be recognized with a period of 4 months and 8 months at I4 (Fig. 18A) and with a period of 4 months at I2 (Fig. 18B). From the cross analysis with the rainfall, we see a high common power during the first semester of year 2 with a period of 4 months and 8 months and with a significant level against red noise with the LOS displacement at I4 and at I2 (Fig. 18). During this episode and this 4-month period, the rainfall is leading the LOS displacement at I4 by about one month (arrow pointing straight down, Fig. 18D) and is in anti-phase with the LOS displacement at I2 (arrow pointing left, Fig. 18E). Interestingly, the XWT reveals that the rainfall is also leading during this episode the LOS displacement near the gas exploitations wells (IZA-M14) but with a very low magnitude compared to the Lussagnet case (Fig. 18F).

LOS displacements in the Izaute zone between December 2015 and April 2016
The times series of the LOS displacement at the monitoring wells I5, I4, I2 and I3 are compared with the rainfall and the Surface Soil Moisture (SSM) during the time interval between December 2015 and April 2016 (Fig. 19). The LOS displacement at I5 is leading the displacement at I4, which is leading the LOS displacement at I2. There is no significant LOS displacement at I3. As indicated by the phase of the XWT (see previous section), the rainfall is leading the LOS displacement at I4 by about one month and the LOS displacement at I2 by about two months. Interestingly, the LOS displacement at I4 displacement is very similar in phase with the LOS displacement at L2, with a leading phase of L2. These results are consistent with the depth of the clay layer which is 3 m at L2, 4 m at I4 and 5 m at I2 given by the geological logs from the BRGM subsurface database (Fig. 20).             where edge effects might distort the picture is shown as a lighter shadow.