Tracking and measuring of clay shrinking and swelling using spaceborne remote sensing

New capabilities for measuring and monitoring are needed to improve the shrink/swell hazard. A new French experimental site at Chaingy (Centre-Val de Loire) has been instrumented using extensometers at the surface and soil moisture sensors in the clay layer. Here we show by direct comparison between remote and in situ data for a period longer than three years that the vertical ground displacements are well-captured by the Multi-Temporal Synthetic Aperture Radar Interferometry (MT-InSAR) technique. In addition to the one-year period, two sub-annual periods that reflect both average ground shrinking and swelling timeframes are unraveled by a wavelet-based analysis. Moreover, we propose here to combine data from two satellite sensors, namely the surface soil moisture from the SMOS satellite with the vertical displacement from Sentinel-1A/1B. The relative phase difference between both time series for the shrinking and swelling periods is calculated for assessing the variations in terms of depth and thickness of the clay layer. An Electrical Resistivity Tomography (ERT) survey is carried out in the studied zone to validate this new method relying fully on remote sensing observations. With regard to future works, the same method coupling different satellites acquisitions may be scaled up to track the expansive clays in a larger area.


Introduction
Among the various natural risks in France, the risk due to shrinking and swelling of subsurface clays is the second most important cause of financial compensation from insurance companies behind the flooding risk. In 2010, the first global shrink/swell hazard map of metropolitan France, based on 1:50,000 geological maps, geotechnical data and spatial distribution of building damages is published by BRGM. So far, in situ monitoring of soil moisture/ground movements and ex situ clay characterization are the traditional way to study this risk, as carried out at the Mormoiron site with a Mediterranean climate (Vincent et al., 2007a;Vincent et al., 2007b). Standard ground displacement monitoring techniques (e.g. extensometers, GNSS) provide information on a very limited number of points within an area. Allowing a higher density of measurement points, the monitoring using Multi-Temporal Synthetic Aperture Radar Interferometry (MT-InSAR) techniques has been intensively developed in the last decades to track land subsidence or uplift related to groundwater extraction or recharge of aquifers around large cities (Declercq et al., 2017;Foumelis et al., 2016;Foumelis, 2012). Conversely, the monitoring of expansive clays using InSAR has been the subject of very few studies until now (Bonì et al., 2018;Burnol et al., 2019;Fryksten and Nilfouroushan, 2019). The major limiting factor to this purpose was the non-availability of relatively high-temporal resolution remote sensing datasets. There is indeed the requirement for fine temporal sampling due to the non-linear behaviour with respect to time of the shrink/swell cyclic behaviour (Raucoules et al., 2009). The launch by the European Space Agency (ESA) of the Copernicus Sentinel-1A/1B satellites enables the systematic data provision, with a 6-day repeat cycle at the equator. Since 2016, a new site characterized by a high shrink/swell hazard level and located at Chaingy, an urban area with a semi-oceanic (i.e. slightly continental) climate, has been instrumented. Two in situ extensometers (EXT1 and EXT2), spaced about 12 m apart, and ten soil moisture sensors at 1.2 m depth are deployed in the studied zone (Fig. 1). Between September 2016 and December 2019, the entire Sentinel-1A/1B archive data from both ascending and descending orbit geometries, 183 and 178 scenes, respectively, were processed using the CNR-IREA P-SBAS (Parallel Small BAseline Subset) service (F. Casu et al., 2014;M. Manunta et al., 2019) implemented on the Geohazards Exploitation Platform (GEP, https://geohazards-tep.eu) (Foumelis et al., 2019). The temporal evolution of ground displacement is compared to the in situ measurements (Fig. 2). We provide evidence that MT-InSAR techniques, in our case the P-SBAS (provided at 90m x 90m spatial resolution), enable accurate measurement of the ground motion cycles due to the clay swelling and shrinking. Moreover, we compare the in-ground soil moistures to the surface soil moisture acquired by the SMOS satellite (Cabot, 2016). Finally, we propose a completely new method for evaluating the depth and the thickness of subsurface clay layers relying fully on remote sensing observations, namely the use of the relative phase difference between Sentinel-1 InSAR displacement and SMOS surface soil moisture time series. We use the electric tomography survey in order to validate this methodology in the studied zone.

Results and Discussion
Intercomparison of InSAR and in situ displacement time series Displacements measured by InSAR are in the direction of the Line-of-Sight (LoS) of the satellite, while displacements measured by extensometers depict motion in the vertical direction. When exposed to moisture, the magnitude of the vertical expansion of soils will depend on the amount of expansive clay minerals in the subsurface. In addition, to the vertical displacement, the subsurface heterogeneities ( i.e. clay depth, thickness and percentage of sand) may introduce spatial variability to the observed ground displacements. We assumed here that Qualitatively, there is a positive correlation between the vertical displacement VERT-ASC-DES at WP as calculated by InSAR using the satellite observations and as measured by groundbased extensometers EXT1 and EXT2 (Fig. 2a). Quantitatively, the least squares correlation coefficient (cc) is much higher for EXT2 (cc=0.69) than for EXT1 (cc=0.47). During the threeyear period, the ground displacement time-series may be decomposed with the addition of a linear trend component T, a seasonal component S and an irregular residual component. There is a linear trend with a shrinking of about -2.2 mm/yr at EXT1 while the measured motion at EXT2 is merely a cyclic component with a shrinking of about -0.3 mm/yr (Fig. 2a). During the same period, the linear trend measured by Sentinel-1A/B lies between both values of EXT1 and EXT2 with a global shrinking of about -0.55 mm/yr . Concerning the cyclic part, the expansion is up to three times higher at EXT2 than at EXT1, with average swelling of 9.4 ± 2.0 mm and 3.1 ± 1.0 mm, respectively (Fig. 2a). Electrical resistive tomography has been carried out in order to image the lithology stratification and the heterogeneity of the clay layer (i.e. depth, thickness, fraction of clay) (see method section). Along three 28.5 m long profiles in the WP cell ( Fig. 1), the resistivity values range from 8 to 100 Ω.m and a three-layer stratified subsurface is highlighted ( Fig. 3 and 4). We used two drill cores along the SC1-SC2 profile to define the clay subsurface material by the resistivity value (Long et al., 2012) (Fig. 4 and Supplementary Fig. 3). A first layer above the depth of 0.59 m ± 0.22 m corresponds to topsoil and backfill with some clay lenses ( Fig. 3 and 4). A second layer with a mean thickness of 1.5 m ± 0.42 m corresponds to the clay material with a resistivity lower than 17 Ω.m ± 1 Ω.m.
Below 2.09 m ± 0.42 m, a third layer corresponds to sandy limestones with a resistivity gradient due to the weathering process. The depth (approx. at 0.5 m) and thickness (about 2.3 m) of the clay layer at EXT1 and EXT2 are very similar (Fig. 3b). However, the resistivity of the clay layer at EXT2 is much lower compared to EXT1, indicating that the clay fraction is much higher beneath EXT2 (Long et al., 2012). That is also consistent with the expansion difference between the extensometers. Over the same period, the swelling magnitude measured by VERT-ASC-DES is 4.5 ± 0.5 mm that is within both extensometer values. This result in agreement with the assumption that the effect of the resolution of the remote sensing method using a 90 m x 90 m cell is an averaging of the vertical displacement in that cell.

Shrinking and swelling periods using Fourier power spectra and continuous wavelet transform
We use first the Fourier spectra of vertical displacements to calculate the main frequency components of time series of both extensometers as well as at WP and East Point (EP) cells of the P-SBAS grid (Fig. 5). The one-year seasonal period is found in the power spectrum of the four displacement times series at EXT1, EXT2, WP and EP (Fig. 6). In the frequency band between 2 and 4 yr -1 (corresponding to a period between 0.5 and 0.25 yr), there are three other frequencies which are noted F1, F21 and F22 (order of increasing frequency) and P1, P21 and P22 periods (order of decreasing period). It should be noted that the frequency values for all displacement time series are close to the first three sub-annual components of the Fourier transform of the soil moistures E1 and W5 at 1.2 m depth (Fig. 6a) and the two ascending and descending time series of the surface soil moisture (SSM-ASC and SSM-DES) as measured by SMOS (Fig. 6b).
We use now the continuous wavelet power spectrum in order to unravel the intermittent physical mechanisms of the shrink/swell process (Fig 6). The continuous wavelet transform (CWT) tool permits the recognition of power in time-frequency space, along with assessing confidence levels against red noise backgrounds (see method section). September 2018 corresponds to the end of the shrinking period for EXT2 and WP cell (Fig. 2a) and there is a maximum of spectrum power at the P1 period (at 2 yr in Fig. 6). December 2017 corresponds to the start of the swelling period for EXT2 and there is a maximum of power at the P21 and P22 periods (at 1.25 yr in Fig. 6). Both swelling periods at WP cell are unraveled at different times, before 1.25 yr for P21 period and after 1.5 yr for P22 period. The same behaviour is observed at EXT1 extensometer and within EP cell (data not shown). In conclusion, the continuous wavelet transform highlights two sub-annual periods that reflect both average ground shrinking (P1) and swelling timeframes (P21).

Intercomparison of in situ soil moistures and SMOS satellite surface soil moistures
We investigated here the soil moisture variations acquired by E1 sensor (near EXT1, see Fig.   1) and W5 (near EXT2) that are both located at the same 1.2 m depth inside the clay layer (Fig.   3). The temporal variations of E1 and W5 are strongly correlated (cc = 0.73) during the threeyear period (Fig. 2b). For the seasonal one-year period, there is a slight leading of E1 moisture relative to W5 (about 0.75 months) as calculated by the cross wavelet transform (XWT) (see method section). Indeed, the infiltration time of the meteoric water inside the clay layer at 1.2 m depth is lower in the E1 case, since the clay fraction of low permeability is shallower at EXT1 compared to EXT2.
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. We investigated also the SSM acquired by the ascending (SSM-ASC) and descending (SSM-DES) orbits of the SMOS satellite, showing rather positive correlation (cc = 0.48). For the one-year period, SSM-ASC and SSM-DES are both in-phase and there is a slight leading of SSM-ASC relative to SSM-DES for the 5-month period (about 0.9 months) ( Supplementary Fig. 2). Hence, the SSM-ASC dataset was considered more appropriate for calculating the phase differences to other time series in the following section. For the one-year period, the phase lead of the surface soil moisture SSM-ASC relative to the in-ground soil moistures E1 and W5 calculated by XWT is about 0.16 yr and 0.22 yr, respectively.
Estimating variations of expansive clay depth and thickness using the time series phase difference The phase difference between the surface soil moisture (SSM) and InSAR displacement time series indicates the time lag between the cause and the effect in the shrink/swell process. We try to verify that depth and thickness of expansive clays are strongly linked to the phase difference between both time series for the shrinkage and swelling periods, respectively. We will use and compare the Fourier and XWT analyses to calculate this phase difference (see method section). The XWT spectra will be high in the time-frequency areas where both CWTs display high values, so this helps identify common time patterns in the two data sets. The XWT permits also the recognition of relative phase in time-frequency space (noted ∆Φ). While P1 and P21 are the shrinking and swelling periods found by the Fourier method (as described above), P1* and P21* are the same similar periods we found using XWT ( Fig. 7 and 8). XWT tool will be used as well to calculate the circular standard deviation of this phase difference (see method section).

No significant difference is found between ∆Φ values of both extensometers EXT1 and EXT2
for the shrinking period P1* and the swelling period P21* (Table 1). This result is consistent with the tomography results which show that both depth and thickness of the clay is indeed comparable at the location of the extensometers (Fig. 3b). For the shrinking period P1*, the phase difference in degrees is three times higher at EXT2 (∆Φ = 107°) than at WP (∆Φ = 34°) (Table 1). This result for the shrinking period is further consistent with the tomographic findings indicating clays lenses in the first 0.5 m subsurface layer of the three 28m-wide profiles in the WP cell ( Fig. 3 and 4). While there is no significant difference for the shrinking period, ∆Φ is higher at EP cell (∆Φ = 110°) in comparison to the WP cell (∆Φ = 72°) for the swelling period (Table 1) For the investigated site, we have shown how coupling SMOS surface soil moisture and Sentinel-1 InSAR displacement has enabled a new understanding of the heterogeneous behavior of the shrink/swell process. The relative phase difference as calculated by XWT with a standard deviation between both time series is used for assessing the variations in terms of depth and thickness of expansive clays. With regard to future works, the same method coupling different satellites acquisitions could be scaled up to track the expansive clays in a larger area and to estimate their relative depth and thickness.

Materials and Methods
Details about the methodologies implemented are given in the Supplementary information file.
Synthetic Aperture Radar (SAR) acquisitions. Copernicus Sentinel-1 mission data were utilized through the Geohazards Exploitation Platform (GEP), which provides access to various Earth Observation satellites repositories (incl. the Copernicus Open Access Hub, PEPS, CREODIAS, ASF and ONDA).
Interferometric SAR (InSAR) processing. The InSAR processing was performed using the Parallel Small Baseline Subset (P-SBAS) algorithm as implemented on the GEP. P-SBAS method is based on the processing of temporal series of co-registered SAR images acquired over the same target area for the generation of Line-of-Sight (LoS) ground displacement time series and velocity maps. The technique allows for the extraction of both linear and non-linear motion components without a priori assumption on the displacement model. Supposing a 1D vertical displacement, the LOS motions are projected to vertical, according to Equations 1-2, whereas combination of different viewing geometries provide us the actual vertical motion component (Equation 3, adapted from (Hanssen, 2001)).
VERT is the vertical displacement, LOS the motion in the Line-Of-Sight direction, θ the incidence angle, and the subscript ASC and DES for the ascending (θASC = 34.4•) and descending track (θDES = 42.8•), respectively.

SMOS Level 3 Surface Soil Moisture (SSM) products.
The SMOS Level 3 SSM products are downloaded through the CATDS dataset site (Cabot, 2016) (https://www.catds.fr). The data are presented over the Equal-Area Scalable Earth (EASE grid 2) (Brodzik et al., 2012) with a sampling of about 25km x 25km and the studied area is included in one grid cell (Fig. 1). We use these 3-day aggregated SMOS-CATDS SSM products for ascending and descending overpasses between September 2016 and December 2019 for each Sentinel-1 acquisition (6 day-repeat cycle).
Signal processing using Fourier analysis. To perform filtering of satellite and in situ signals, time series of measurements are smoothed using a one-dimensional convolution approach with a HANNING window. To carry out the spectral analysis, the filtered time series is first padded with trailing zeros to a length of 100 yr before computing the discrete Fourier transform (DFT) using the fast Fourier transform (FFT) Matlab function. The frequency maxima of Fourier power spectrum are therefore computed with a precision of 1/100 yr -1 . The magnitude and the phase angles of complex FFT values are calculated by Matlab ABS and ANGLE operators, respectively. Whenever the jump between consecutive angles is greater than or equal to π radians, UNWRAP function shifts the angles by adding multiples of ±2π until the jump is less than π.
Signal processing using wavelet analysis. The wavelet transform is especially suited to extract features from low signal-to-noise ratio time-series and to identify localized intermittent periodicities. The Morlet wavelet is adapted to geophysical time series as described by (Torrence and Compo, 1998) who developed the software provided at http://paos.colorado.edu/research/wavelets. We used the Matlab version as adapted by (Grinsted et al., 2004) and provided at http://www.glaciology.net/waveletcoherence. The continuous wavelet power spectrum CWT expands time-series records into time/frequency space. The time-series input data must be equally spaced in time. Although Copernicus satellites have a regular revisit interval (6 days for Sentinel-1A/B and 3 days for SMOS), some acquisitions may be missing or excluded from processing. These missing values are linearly interpolated using a constant time interval. The other time-series data of extensometers and soil moistures sensors are down-sampled using the same time interval. Two individual CWTs can be combined by using the cross wavelet transform (XWT) tool which is computed by multiplying the CWT of one time-series by the complex conjugate of the CWT of the second time-series. XWT image is the 2-D representation of the absolute value and the phase of the complex number in the time-frequency space. For many geophysical phenomena, an appropriate background spectrum is red noise (increasing power with decreasing frequency). ANGLEMEAN function calculates the mean angle found by XWT during a time period and the sigma which corresponds to a standard deviation.
Clay layer characterization using electric method. The data were collected using the IRIS Syscal Pro Plus multi-electrode imaging system with internal multiplexer. Three profiles were acquired using 96 stainless electrodes spaced every 30 cm. Length of the profile was 28.5 m. The measurement was carried out using Dipole-Dipole and Wenner configuration Due to good coupling, no data was removed during pre-processing. Inversion will be carried out using Res2DINV software with a L1 regulation norm. Root means square values was lower than 1% after fives iterations for the three profiles.       In the bottom, the phase lags of SSM-ASC relative to the same times series are shown for the shrinking (P1*) and swelling (P21*, P22*) periods found using the thick contours of XWT that designate the 5 % significant level against red noise.