Wavelet-based analysis of ground deformation coupling satellite acquisitions (Sentinel-1, SMOS) and data 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 double gas storage site (Lussagnet and Izaute) in Southwestern France. Despite quite low vertical displacements (between 4 and 8 mm) compared to the noise level, the cyclic motion reflects the seasonal variations due to charge and discharge during summer and winter periods, respectively. We can simulate the ground deformation at both storage sites by a simple mechanical model. However, ground movements of low-magnitude may be also induced by natural factors, such as the temperature or the soil moisture. Using a wavelet-based analysis, we show there is a soil expansion in the Lussagnet zone that contrasts both in phase and period with the seasonal deformation and that is linked to the surface soil moisture 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 reveals the combination of two different processes driving the ground displacement with the same order of magnitude (about 6 mm), namely the pressure variation of a deep gas reservoir and the swelling/shrinking of the shallow subsurface.

The double gas storage reservoir isobaths between the river L'Adour at the West and La Midouze at the East are shown using red lines and the ~25 km SMOS cell is shown by a black square. On the western edge of Le Houga city is located the Lussagnet well L1 (red circle) and on the eastern edge the Izaute well I1 (red circle). The Latrille piezometric borehole and the Mont-de-Marsan weather station are located at a distance of about 3 km and 27 km from L1 well, respectively.

II/ Sentinel-1 satellite acquisitions using SqueeSAR ® technique
Methods SqueeSAR ® is a second-generation PSInSAR algorithm, a method based on the processing of temporal series of co-registered SAR images acquired over the same target area. The main idea behind PSInSAR is to identify point-like targets (corresponding to single pixels or groups of a few pixels) exhibiting good phase coherence over the entire observation period via proper statistical analyses. Unlike the PSInSAR approach, the SqueeSAR ® technique allows the measurement of surface displacements by exploiting both point-wise coherent scatterers (i.e. the PS) and partially coherent Distributed Scatterers (i.e. the DS).
The sensor to target direction is orthogonal to the orbit and it is called the Line of Sight (LOS).
As any other PInSAR technique, SqueeSAR ® measures the projection of the real displacement along the LOS, that is inclined with respect to the vertical axes. The main implications of the LOS measurements are the following: • 115 images are processed during this 3-year period and 9 images are missing or excluded from processing (see Table S1 and Fig. S2).
The trend part T is approximated by a second order polynomial function of time and the seasonal part S itself may be split into one bottomhole pressure component P and one surface temperature component ST: Our additive decomposition model estimates therefore five coefficients, i.e. both coefficients ( , ) of seasonal part and three coefficients of the quadratic trend part, by minimizing the residual sum of squares (RSS) for all measurement points. IV/ Ground deformation modeling during each gas exploitation cycle

Methods
The vertical displacement ( ) on the ground surface is expressed by an integral equation over the volume where pore pressure changes. For a uniform distribution of pressure change (Δ ) within a horizontally laid disk-shape reservoir (radius of 0 and a depth between ℎ and ℎ ) in a homogeneous, semi-infinite elastic medium, we have the integral equation given by Geertsma (1973): where cylindrical coordinate ( , ) is taken (Fig. S4a), and 0 ( ) is the Bessel function of the first kind for integer order 0. The horizontal component ( , ) is written by a similar integral equation as well. 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 (Fig. S4b). The assumed model parameters are summarized in Table S2. This is a phenomenological, namely kinematic, expression without solving hydrogeological or geomechanical system of the reservoir. From the geological map, we estimate the extension of the reservoir as 5.4 km 2 for Lussagnet and 13.1 km 2 for Izaute, respectively. 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 (SE3) concentrates around the reservoir surface ( Fig. S4c-e). From this estimation, the vertical displacement on the ground surface arrives at 8 mm and 4 mm in the Lussagnet case and in the Izaute case, respectively (Fig. S4c). The deformation pattern is quite independent, meaning that the deformation amplitude is close to zero half way in-between the two reservoirs. As it is clear from Equation (SE3), the amplitude of deformation depends strongly on the material constants ( , , ) given in Table S2. 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 (see main text, 6.4 ± 0.8 mm). It can also be noted that a smaller reservoir geometry shall lead to a larger amplitude at the center of the reservoir.
The root mean square (RMS) of LOS displacement in the reference zone is 0.94 mm (Fig. 2a). Supposing a sinusoid shape, the peak-to-peak LOS displacement is 2√2 × 0.94 = 2.66 mm. This value is our estimated noise level in the DInSAR time series in the studied zone.
Around the gas exploitation wells (L1 and I1), the predicted horizontal displacement is less than 2 mm (Fig. S4d-e). This value is less than the estimated noise level in the reference zone (2.7 mm). Therefore, we can neglect the contribution of horizontal displacement around L1 and I1 wells and transform the displacement from LOS to the vertical using = cos ⁄ .
Around the monitoring wells (L2-4 and I2-5), the predicted vertical displacement during each cycle is less than 2.5 mm (Fig. S4c). This value is less than the estimated noise level in the reference zone (2.7 mm). Therefore, we investigate by a wavelet-based analysis the LOS displacement at these monitoring wells to track another process driving the ground deformation in the studied zone.  V/ Surface Soil Moisture data extracted from SMOS level 3 products

SMOS satellite was successfully launched on the 2n d of November 2009 by the European Space
Agency (ESA). SMOS uses microwave radiometry for estimating soil moisture. We use here the term Surface Soil moisture (SSM) to refer to the volumetric soil moisture in the first few centimeters (0-5 cm) of the soil. L-band radiometry is achieved resulting in a ground resolution of 50 km. SMOS Levels L0 to L2 processing products are developed by the ESA. The products are divided in half orbits, from pole to pole, ascending or descending, spanning about 50 minutes of acquisition 1 . Level 3 products are geophysical variables with improved characteristics through temporal resampling or 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) 2 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. A recent published evaluation of the accuracy of SMOS-CATDS products between 1 January 2016 and 30 June 2017 in southwestern France 3,4 gives an average bias, RMSD (Root Mean Square Difference), ubRMSD (unbiased Root Mean Square Difference) of about -10 vol.%, 12 vol.% and 7 vol.%, respectively. The SMOS-CATDS products thus underestimate the surface soil moisture (with quite high negative bias values) but the ubRMSD of 7% is similar to the best values obtained by other satellites in Southwestern France (about 6%) 3,4 . The filtered daily global maps are aggregated over different time periods. The SMOS revisit period is 3 days at the equator, so the Earth's surface is covered by SMOS field of view in 3 days and a 3-day global product with a 3-day moving window is provided by CATDS by aggregating the daily global maps. The retrieved SSM data are filtered using the Data Quality Index (DQX) provided in the product (DQX < 0.06 m3/m3). The filtering processor is then used over the Level 3 retrievals, to select the best estimation of soil moisture if several retrievals are available for a given day.

Results
From the analysis of the CWT of the rainfall time series (Fig. S5), we can observe two cycles (8-month and 4-month) during the year 2 in addition to the 1-year cycle. The CWT of the Surface Soil Moisture is similar to the rainfall CWT for the 4-month period, with a higher power for the 1-year cycle but with less power for the 8-month period. From the CWT analysis of the Latrille piezometric level, we can observe only the 1-year cycle and no patterns in higher frequencies.